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Abstract 


A stochastic-geometric landsurface reflectance model is formulated and tested 
for the parameterization of spatially variable vegetation and soil at subpixel scales 
using satellite multispectral images without ground truth. Landscapes are 
conceptualized as three-dimensional Lambertian reflecting surfaces consisting of 
plant canopies, represented by solid geometric figures, superposed on a flat soil 
background. Multiple scattering among landsurface components is neglected. The 
model is cast within the framework of an existing theoretical model of upwelling 
solar radiance for optically-thin atmospheres, as observed by a nadir-viewing 
satellite. 

A computer simulation program is developed in order to investigate image 
characteristics at various spatial aggregations representative of satellite 
observational scales, or pixels. In particular, the evolution of the shape and 
structure of the red— infrared space, or scattergram, of typical semivegetated scenes 
is investigated by sequentially introducing model variables into the simulation. The 
correlation between canopy and shadow is identified as a principal mechanism 
contributing to the frequently observed tasseled cap of red— infrared scattergrams of 
semivegetated landscapes. A Sampling Scale Ratio is formulated as a quantitative 
criterion that identifies when that correlation occurs. 

The analytical moments of the total pixel reflectance, including the mean, 
variance, spatial covariance, and cross-spectral covariance, are derived in terms of 
the moments of the individual fractional cover and reflectance components. The 
moments are applied to the solution of the inverse problem: The estimation of 
subpixel landscape properties on a pixel— by— pixel basis, given only one 
multispectral image and limited assumptions on the structure of the landscape. The 
inverse procedure involves the formulation of conditional moments for subsets of 
pixels that possess similar properties, and that can be identified through their 
common orientation in red-infrared scattergrams. The analysis is facilitated by 
assuming geometric similarity among canopy elements and by assuming a functional 
relationship between fractional covers in the case of large Sampling Scale Ratios. 

The landsurface reflectance model and inversion technique are tested using 
actual aerial radiometric data collected over regularly spaced pecan trees, and using 
both aerial and Landsat Thematic Mapper data obtained over discontinuous, 
randomly spaced conifer canopies in a natural forested watershed. For the Landsat 
case, adjacency effects are neglected by assuming low interpixel contrast. Different 
amounts of solar backscattered diffuse radiation are assumed and the sensitivity of 
the estimated landsurface parameters to those amounts is examined. 
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Chapter 1 
INTRODUCTION 



The purpose of this research is to formulate and test a physically-based 
reflectance model that characterizes the spatial variability of multispectral images 
obtained over semivegetated landscapes. The immediate objective is to develop a 
flexible, physically-based algorithm for estimating the amount of subpixel 
vegetation cover (i.e. horizontal fractional cover) on a pixel-by-pixel basis using 
only one set of satellite multispectral data under clear-sky conditions, without 
ground truth, and without having to compute the numerous scattering and 
absorption parameters that govern atmospheric radiative transfer. The focus is on 
natural landscapes that exhibit random behavior in the size and location of 
individual plants, and in the soil background reflectance. The goal is to 
accomodate both the subpixel scales associated with the bulk physical properties of 
the plant canopy (overall geometry, height, and diameter) as well as the regional 
scales associated with the variability in soil background reflectance. The 
long-term objective is to provide a framework for the physically-based 
parameterization of mesoscale landsurface hydrology using remote multispectral 
observations. 



The physically-based parameterization of the large-scale coupled heat and 
moisture fluxes of semivegetated landscapes is a unsolved problem in hydrology 
(NASA, 1988). The principal difficulty arises from the complex spatial and 
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temporal heterogeneity of the interrelated hydrological, geophysical, and biological 
processes, as well as the many landsurface variables which define them. The 
hydrologic variables, including soil moisture, vegetation type and amount, 
hydraulic conductivity, and temperature, often exhibit random behavior and 
possess spatial scales of variation that are much smaller than the overall scale of 
the parameterization (mesoscale or greater). Such heterogeneity limits the fidelity 
of the usual homogeneous mathematical description of individual hydrologic 
processes (evapotranspiration, infiltration and runoff) and engenders the formidable 
logistical problem of how to acquire regionally representative hydrologic data in an 
efficient, cost-effective manner. 

Most classical approaches to mesoscale investigations, such as flood 
forecasting, river basin management, and environmental impact assessment, rely on 
standard hydrometeorological data obtained at ground-based stations. Such 
stations, usually few in number, are often located at airports or agricultural 
research sites and almost never in remote watersheds. Regional surface fluxes are 
generally estimated using an area-weighted extension of point estimates calculated 
for the ground station data. However, since the ground stations do not necessarily 
capture all of the basinwide spatial variability in soil, vegetation and climate, and 
since the location of a station itself may be biased toward a particular hydrologic 
regime, the accuracy of that approach is highly uncertain. 

The problem is more critical in global— scale hydrologic parameterizations, as 
modeled within atmospheric general circulation models (GCMs). Those models 
possess typical grid scales (100 km) which are much greater than the spatial scale 
of fluctuation of the hydrologic processes themselves (i.e. 1-100 km, See 
Smagorinsky, 1978). Although, traditionally, many GCMs prescribed simple 
surface flux models (i.e. Manabe, 1969; Washington and Parkinson, 1986), more 
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detailed algorithms have been recently proposed which include a large number of 
soil and vegetation properties (e.g. Dickinson, 1984; Rind, 1984; and Sellers et al, 
1986). However, the practical benefit of those newer, more sophisticated 
parameterizations is not realized due to the lack of spatially detailed global 
landsurface hydrologic data. 

A specific need for more detailed spatial landsurface data may be beneficial 
for estimating subsurface properties, through the application of a time-averaged 
one-dimensional statistical-dynamic representation of the climate-soil-vegetation 
system. Eagleson (1982) mathematically formulated three interrelated hypotheses 
describing the short, medium, and long-term equilibrium states of soil, vegetation, 
and climate for natural undisturbed systems. He argued that in water-limited 
systems, the short-term ecological pressure minimizes water demand stress through 
adjustment of vegetation canopy amount, the medium term pressure selects plant 
species for minimum water use, and in the long term, vegetation and climate 
modify soil properties in the root zone in a synergistic manner to reach a 
climatic-climax state such that biomass productivity is maximized. Limited 
testing of those hypotheses have been made on several catchments (Eagleson and 
Tellers, 1982) and on two savanna systems (Eagleson and Segarra, 1985). 

Estimation of Soil Hydraulic Properties. Jasinski (1987) proposed to apply 
the equilibrium hypotheses to estimate soil hydraulic properties of natural 
water-limited systems, using vegetation density estimated from satellite data. A 
critical step in this analysis was the development of an algorithm for estimating 
spatially-variable fractional canopy cover, which controls the bounds of the 
rejection probability of the hypothesis. He collected extensive data for that 
investigation which is currently underway (See Jasinski and Eagleson, 1986; 
Eagleson and Jasinski, 1988). 
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1.3 AppliCfttiQn of Sa tgllitq Radiometric Observations to Regional Hydrology 


The critical need for detailed spatial and temporal knowledge of regional 
landsurface processes has led researchers to investigate the application of 
electromagnetic radiation data obtained from satellite platforms. Those data 
consist of instantaneous, spatially integrated observations of the electromagnetic 
radiation fluxes emitted or reflected from the earth's surface and atmosphere. The 
reasonableness of this approach rests not only in the relative facility of covering 
extensive areas in a matter of seconds, but also in the fact that many of the 
physical properties which describe the heat and moisture fluxes (i.e. landsurface 
and atmosphere composition, temperature) also govern radiative tr ans fer 
Unknown landsurface and atmospheric properties are theoretically determined by 
solving the inverse problem. That is, given a set of remote electromagnetic 
observations, one inverts the radiative transfer equation using a particular 
wavelength, or combination of wavelengths, so that the parameters of the 
reflecting surface and the medium can be retrieved. The challenge lies in carefully 
choosing specific wavelengths within the spectrum which respond to the presence 
(or absence) of the particular object or constituent under investigation. 

Although the theoretical radiative transfer aspects of remote sensing are well 
understood (i.e. Chandrasekhar, 1960; American Society of Photogrammetry 

Slater, 1980; Tsang et al, 1985), because of the variability and large 
number of landsurface and atmospheric parameters, and the limited number of 
satellite systems, the use of remotely sensed data alone has not been sufficient to 
estimate landsurface fluxes as of date. Successful retrieval of individual 
landsurface parameters has been achieved (See ASP, 1983), however, which can be 
incorporated into classical surface heat and moisture flux equations. For example, 
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the microwave portion of the spectrum has been shown to be sensitive to soil 
moisture content because of the large contrast between the dielectric properties of 
liquid water and dry soil (Schmugge, 1983). The visible and near-infrared 
portions of the spectrum are sensitive to vegetation amount due to the contrast in 
reflectance behavior between soil and healthy green vegetation (Colwell, 1974). 
Thermal infrared wavelengths have been shown to be sensitive to surface and 
atmospheric temperature (Chahine, 1983). The published results of those 
investigations are too numerous to mention, although an excellent summary of the 
response of different wavelengths to particular landsurface properties is provided in 
ASP (1983). 

Most evapotranspiration investigations which have been reported typically 
combine one or more remotely sensed landsurface variables, such as albedo and 
temperature, with ground based meteorological data, and apply them to classical 
sensible and latent heat flux terms in the energy balance equation. Several studies 
over relatively homogenous agricultural areas have employed either airborne or 
hand-held radiometers. Camillo et al (1983) developed an energy and moisture 
balance model of the upper soil and lower atmosphere for use with remotely sensed 
surface temperature and soil moisture and standard meteorological data, and 
applied the model over bare soil. A similar approach was later used by Gurney 
and Camillo (1984) over wheat and barley. Van de Griend and van Boxel (1989) 
extended the approach to include multilayer canopy representation, and tested the 
model over maize. Reginato et al (1985) combined remotely sensed reflected solar 
radiation and surface temperatures using airborne sensors, with ground 
meteorological data to calculate net radiation and sensible heat flux over wheat. 

Theoretically, the same approach as above can be used for investigating 
surface fluxes using satellite-based sensors. However, practical application has 
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been limited due to complexities arising from atmospheric effects and the problem 
of characterizing landsurface heterogeneity, especially with regard to vegetation 
and soil moisture. Most hydrologic studies using remote sensing data limit the 
study area to spatially homogeneous landscapes, or they assume spatial 
homogeneity. Price (1982) applied visible and thermal infrared data from the 
Heat Capacity Mapping Mission (HCMM), to the energy balance equation to 
estimate landsurface thermal inertia and moisture availability over principally 
grassland areas in Washington and Oregon. Similar approaches have been used by 
Carlson (1986) using GOES data over grasslands in Kansas, and by Taconet et al 
(1986) using AVHRR data over wheat, and by Pierce and Congalton (1988) using 
simulated TM data over mixed conifer forests in California. 

While it is evident from above that significant progress has been achieved in 
the analytical treatment of satellite data for regional hydrologic investigations, 
there is substantial need for additional research. Ongoing programs within the 
International Satellite Land-Surface Climatology Project (ISLSCP), such as the 
First ISLSCP Experiment (FIFE) (Sellers et al, 1988; Hall et al 1989) and the 
Hydrological Atmosphere Pilot Experiments (HAPEX) (Andre, 1986) are currently 
bringing into focus the current limitations, and research needs, for applying 
satellite data to the study of landsurface hydrological processes. Research areas 
include developing methods to improve extraction of all the information hidden in 
the multispectral data, as well as rethinking the basic structure of classical energy 
and moisture balance equations to accommodate the particular characteristics of 
satellite data. 

A recurrent problem identified by many of the above studies is the 
characterization of the spatial variability of vegetation amount, whether it be 
fractional cover, biomass, or leaf area. Knowledge of vegetation cover is 
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hydrolog i cally important with regard to the partitioning of bare soil evaporation 
and transpiration, as with its influence on albedo, moisture storage, and surface 
temperature. Currently, however, no physically-based technique exists for 
estimating vegetation amount for large areas using only multispectral data. Most 
estimation procedures are empirical, involving ratios, differences, or 
transformations of signals in red and near infrared bands (Perry and 
Lautenschlager, 1984). Such indices respond to changes in vegetation amount 
primarily due to i) the relatively high radiation absorption capacity of chlorophyll 
in the red band, and ii) the high reflectance properties of the leaf structure in the 
near infrared band, as compared to soil. Despite that sensitivity to vegetation 
amount, however, vegetation indices provide only limited understanding of the 
physical structure of the scene. They generally require a large number of training 
samples and can exhibit inordinate scatter for equivalent amounts of vegetation. 


Estimation of Subpixel 


The physically-based characterization of spatially variable plant cover, using 
satellite multispectral data, is a critical constraint to hydrologic parameterization 
in many semivegetated landscapes. An important example is the natural semiarid 
region of most of the southwest United States, which typically consists of a 
random distribution of different size tree or shrub canopies interspersed with a 
mixture of grasses and bare soil. Also important are agricultural lands during 
early growth stage, in which the crops are more uniform in size and spacing, and 
separated primarily by bare soil. 

One of the major problems in trying to use satellite multispectral data in 
semivegetated regions is that the plant canopy or size typically varies at 


24 


characteristic scales (several meters) much smaller than the spatial resolution of 
current satellite pixels (several tens of meters). At the same time, soil 
background reflectance varies over a wide range of length scales (meters to several 
thousand meters) due to geoclimatic factors affecting its physical structure and 
chemical composition. Since satellite observations integrate the reflectance of all 
elements within the pixel, subpixel information such as fractional cover, leaf area, 
surface roughness, is apparently lost. Thus, techniques are needed to disaggregate 
the individual subpixel components, while accounting for the variations in both 
soil and vegetation reflectance. 

The above problem of estimating subpixel vegetation cover can be 
graphically illustrated using Figures 1.1 through 1.4. Figure 1.1 depicts a 
hypothetical landscape viewed through a red filter, in which the plants appear as 
dark circular disks of constant reflectance and of random size and distribution. 

The plants are superimposed on various soils which can possess different 
reflectances as indicated by the different shading on the figure. Also drawn on the 
scene is the outline of nine satellite pixels, each of which contains any number of 
trees or any type of soil background. 

Figure 1.2 illustrates the remotely-sensed image of the scene in Figure 1.1 
using the same red filter. Only one reflectance value exists for each pixel, 
obtained through the sensor's spatial integration of the reflectances of all the 
elements in the pixel. Thus, the subpixel information (size, number and 
distribution of the trees, distribution of soil reflectance) can not be discerned. 

The same scene as viewed in the near infrared spectrum is depicted in Figure 
1.3. In that band the plants generally appear brighter than the soil. The soil 
reflectance, although brighter than in the red region, possesses approximately the 
same spatial distribution. The remotely-sensed image of the near infrared scene 
also averages out the components of the pixel, as shown in Figure 1.4. 
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Figure 1.1 Hypothetical semivegetated scene Figure 1.2 Remotely-sensed image of 

viewed through red filter. scene in Figure 1.1. 
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Figure 1 .3 Hypothetical semivegetated scene Figure 1 .4 Remotely-sensed image of 

viewed through near-infrared filter. scene in Figure 1 .3. 





The inverse problem can be summarized in the following questions. Given 
only the two remotely-sensed images, to what extent can the subpixel properties 
of the original scene be retrieved? Can a physically-based inverse method be 
developed that exploits the true multispectral nature of the data (i.e. two 
observations per pixel), by incorporating cross spectral correlations? Are there 
assumptions with regard to the geometry or spatial distribution of the plants and 
the soil, which may be useful toward the solution? Also, are there particular 
relationships between the soil and vegetation which preclude certain solutions, such 
as the relationship between soil reflectance and vegetation amount, or vegetation 
amount and shadow? Those questions need to be considered to adequately address 
this problem. 

In the real world, the problem is still more complex than the preceding 
example, especially for large-scale parameterizations. In addition to subpixel 
variations, one encounters many regional-scale variations in soil and vegetation 
reflectance due to a variety of geoclimatic factors. For instance, changes in slope 
and aspect induce corresponding changes in scene reflectance through an effective 
altering of the illumination and viewing angles. Changes in elevation, slope, and 
aspect also cause scene variability through their indirect influence on such 
properties as soil moisture, and vegetation species and density. Vegetation 
reflectance can change with plant size and density, and with changes in underlying 
soil reflectance. In such regions the stochastic nature of the vegetation and soil 
properties must be accommodated. Thus, the problem consists of trying to 
discern, through the interpretation of multispectral data, not only the small-scale 
(i.e. subpixel) variability, but also the regional-scale correlations that might exist. 

The ability to estimate regionally variable subpixel vegetation cover would 
be invaluable to the solution of the large-scale parameterization problem. 
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Successful development of such an algorithm has potential utility for investigating 
not only vegetation amount but also other surface geophysical properties requiring 
subpixel resolution. 

1.5 Physical Basis of Red-Infrared Scattergrams: Preliminary Analysis 

A common attribute of the red-infrared data spaces of vegetated images is 
that they often take on a triangular shape, or tasseled cap. The physical basis of 
this shape was first investigated by Kauth and Thomas (1976) using the Suits 
model (1972) applied to a homogeneous layer of crops. They explained the 
seasonal progression of the tasseled cap in terms of the growth, maturation, and 
senescence of crops. 

It was shown during the initial stages of this research that individual images 
of many natural semivegetated landscapes also exhibit triangular shapes when 
plotted in the red-infrared space (Jasinski, 1987). Three such scattergrams, 
constructed using three segments of Bands 2 and 4 of a Landsat 2 Multispectral 
Scanner (MSS) image in the vicinity of Taos, New Mexico, are shown in Figures 
1.5 through 1.7. The first triangular scattergram, Figure 1.5, consists of a plot of 
the red-infrared data pairs of all the pixels in a segment covering about 400 
square kilometers. The region contains a variety of semivegetated landscapes, 
ranging from bare soils, grasses shrubs along the flat lands in the valley, to 
pinyon-juniper in the foothills, to ponderosa pine and douglas firs in the 
mountains. The resulting characteristic triangular plot is generally flat at the 
base, but curved along the top. Figure 1.6 contains a smaller segment of pixels 
covering about 25 square kilometers, located principally in the foothills region 
dominated by pinyon-juniper trees. The resulting shape is also triangular in form, 
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Figure 1.5 Plot of MSS Band 4 versus Band 2: Entire Taos Study Area 
(40,000 pixels). 
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but much smaller than in Figure 1.5 and with data pairs exhibiting. overall lower 
values. Finally, Figure 1.7 possesses pixels obtained primarily over 80 square 
kilometers of grasslands. The scattergram again exhibits a triangular shape, but 
with a wider base than from the foothills region. 

Although triangular in shape, the above scattergrams possess many features 
not present in the Kauth-Thomas investigation. For example, they were 
constructed from natural vegetated regions and from only one MSS image. 

Further, individual pixels of the scene possessed different amounts of vegetation 
species, cover, soil and topography. Thus, the physical mechanisms which lead to 
those triangular shapes can not be explained simply in terms of the season growth 
of plants as noted by Kauth and Thomas. 

Using a simple area-weighted combination of vegetation and soil reflectances, 
it was demonstrated that the characteristic triangular shapes of red-infrared 
scattergrams could be explained in terms of percent canopy cover, variable soil 
and vegetation reflectance, and shadows (Jasinski, 1987). For instance, using 
typical values of vegetation and soil reflectance, the hypothetical data space of a 
landscape possessing constant vegetation reflectance, variable soil reflectance, and 
variable canopy cover is, in fact, a triangle (see Figure Al.l). The effect of 
changes in either the vegetation or soil reflectances causes an equal response in the 
shape or position of the triangle (Figures A1.2 and Al.3). The inclusion of 
shadows cast by conical figures causes the triangle to take on a tasseled cap 
(Figure A1.4). By assuming a linear relation between vegetation reflectance and 
percent cover, the data space takes on curved sides (Figure A 1.5). Other 
variations in plant or soil parameters, such as leaf area index and soil moisture, 
can be shown to cause similar effects on the red-infrared scattergram. 

The above exercise demonstrated that a key to solving the inverse problem 
lay in an understanding of the physical basis of the red-infrared scattergram. It 
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further showed that the problem of estimating subpixel vegetation cover was 
inherently tied to the distribution of other landscape properties as well, in 
particular, soil background reflectance and shadows. That early discovery 
provided impetus to the direction of the methodology followed in this thesis. 

1.6 Elaboration of Goals and Methodology 

The research in this report addresses the parameterization of vegetation and 
soil of natural semi-arid regions using multispectral data, for application in 
large-scale hydrology. Specific goals and methodology are summarized below. 

The first step is the formulation of a physically-based canopy-soil reflectance 
model that characterizes the spatial variability of multispectral images of 
semi vegetated landscapes at both subpixel and regional scales. A 
stochastic-geometric reflectance model is developed for that purpose, as it is 
flexible enough to represent the possible stochastic nature of the distribution of 
vegetation and soil, but not so detailed as to render its practical application 
infeasible. The landsurface reflectance model is cast within the framework of an 
existing theoretical model of upwelling radiance observed by a nadir-viewing 
satellite. 

The second step is the development of an understanding of the shape and 
structure of red-infrared scattergrams in terms of the physical properties of the 
landscape. That is achieved through computer simulations of idealized 
semi vegetated landscapes using the above canopy-soil reflectance model. The 
principal mechanisms that contribute to the evolution of the triangular shape of 
red-infrared scattergrams are identified through the sequential altering of a given 
variable in the model. 
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Next, an inverse method is developed for estimating fractional (i.e. subpixel) 
canopy cover on a pixel-by-pixel basis using only multispectral data without 
ground truth. The inverse procedure chosen for the analysis is the method of 
moments in which the theoretical moments of the canopy-soil reflectance model 
are derived and equated to the actual moments of the multispectral image. The 
procedure makes use of the understanding of the shape and structure of the 
red-infrared scattergram, as it requires writing conditional moments for portions of 
the scattergram where the pixels are assumed to possess one or more similar 
attributes. The inverse procedure is further facilitated by assuming geometric 
similarity among the plant and shadow geometry, and through the formulation of 
a functional relationship among the fractional covers based on the similarity 
parameter and an assumed spatial distribution of the plants. A Sampling Scale 
Ratio is developed, based on the relative scales of the cover components relative to 
the scale of the pixel, which quantifies when the correlation among fractional 
covers exists. 

After testing the inverse method on the simulated images to determine its 
effectiveness, the inverse procedure is then tested on actual aerial and satellite 
multispectral data. The field sites include one agricultural area and one natural 
forested scene, both in Arizona. Other properties of the landscape, such as the 
spatial distribution of soil reflectance, are also examined at that time. 
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Chapter 2 


RADIANCE OBSERVED BY SATELLITE 


This chapter briefly reviews existing fundamental radiative transfer theory 
related to remote sensing of landsurface properties from a satellite platform. The 
purpose is to provide an analytical framework, using an existing model of 
surface-reflected and atmospheric-scattered radiance, into which the canopy-soil 
reflectance model developed in Chapter 4 can be incorporated. While the 
principles in the present chapter generally apply to all wavelengths and media, the 
emphasis is on visible and near-infrared wavelengths interacting in the earth's 
atmosphere. Particular characteristics relevant to Landsat Multispectral Scanner 
(MSS) and Thematic Mapper (TM) sensors are reviewed. 

2.1 Radiative Transfer Theory 
2.1.1 Definitions and Nomenclature 

The propagation of electromagnetic energy of a particular wavelength, A, is 
described in terms of the specular radiance, L (A), defined as the specular radiant 

3 

energy flux per unit wavelength per unit solid angle per unit of projected area in a 
specified direction. A graphical illustration of radiance emerging from an 
elemental area, dA, centered at P(x,y,z) is drawn in Figure 2.1. The sketch 
depicts a portion of specular radiant flux, d({) (A), that emerges from dA and that 
is confined within the solid angle, dQ. The flux makes an angle 9 with the 
vertical z axis, and an angle ip with the x axis in a clockwise manner when viewed 
from the bottom. 
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Mathematically, the relation between specular radiance and radiant flux is 
given by, 


L (A) = d2( ft( A l 
' dft d A c o s 0 


( 2 . 1 ) 


where the units of radiant flux are Watts per micrometer (W /xm -1 ) and those of 
radiance are Watts per square meter per steradian per micrometer (W m - ^ 

sr _ V m_ *)- 

Integration of the specular radiance over a finite bandwidth yields the total 
radiance within that interval of the spectrum, or 


L (A) = J A ^L s (A')dA' 


( 2 . 2 ) 


where now A represents the center of the band, and the units are Watts per 

—2 —1 

square meter per steradian (W m sr ). While this measure is also specular in 
the sense that it represents the radiation only within a particular bandwidth, for 
simplicity, the adjective "specular" is deleted as the wavelength dependency is 
adequately represented by the nomenclature. 

The integration of the radiance emitted from dA over the entire spherical 
angle is termed radiant exitance, M(A), or 


M(A) = f L(A) cos0 dft 
2ft 


( 2 . 3 ) 


It possesses units of Watts per square meter per micrometer (W m ^/zm *). 
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The total radiance impinging on a surface from all directions can be defined in an 
analogous manner and is termed irradiance, E(A), having the same units. 

A beam of radiation impinging on a surface at a given angle can be 
reflected, absorbed, and transmitted. The relative amounts of each component, as 
well as the direction of the reflected radiation, depend on the properties of both 
the incident beam (i.e. wavelength, angle of incidence) and the properties of the 
surface (i.e. roughness and chemical composition, etc.). Figure 2.2 illustrates the 
geometry describing reflectance due to a source beam, L'(A;0\^') impinging on dA, 
as viewed from an observer located at (0,tp). The observed beam, dL(A ;9,ip) 
represents only a portion of the total radiation reflected from the surface, as 
radiation is simultaneously being scattered in other directions. This angular 
dependent reflectance, due to a single source, is termed bidirectional reflectance, 

It is a dimensionless quantity and is mathematically defined, 


dR(A 


dL(A:0.t/>)cosfl dfl 
V(A;>,^')cos0'dfr 


(2.4) 


When the incident radiation arrives from sources throughout the entire 
hemispherical angle, a useful quantity is the bidirectional reflectance distribution 
function (BRDF), or 


f(A;<M';M) 


Ww7V) 


(2.5) 


with units of inverse steradians (sr - *) (Slater, 1980). 

When the reflected radiation is independent of angle, the surface is described 
as Lambertian. In many practical remote sensing investigations over natural 
landscapes, surface reflectance can often be assumed as Lambertian, especially at 
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Source beam 


Observed beam 



Figure 2.2 Geometry of surface reflectance 
(After Slater, 1980). 
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low zenith angles (Slater, 1980). One exception is water which exhibits 
bidirectional reflectance. The specular reflectance, R(A), of a Lambertian surface 
is simply, 


“ ETaT* (2-6) 

(Slater, 1980). It is also a dimensionless quantity. 

2.1.2 Radiative Transfer Equation 

A beam of radiation traversing a medium can be attenuated due to 
absorption and scattering, or enhanced by emission and multiple scattering, as it 
interacts with the medium. A graphical sketch of this process is shown in Figure 
2.3. For atmospheric applications, it is often sufficient to assume a plane-parallel 
medium in which the stratification lies along the z axis (vertical). In such cases 
the radiative transfer equation in its most general form is, 

» I UXrjiA a = L(A, - J(A,t- A *P) (2.7) 

where J represents the source function (W nf^sr^/mf 1 ), fi equals cos 0, and r is 
the optical thickness (dimensionless). The source function characterizes all the 
scattering and emission from external sources which enhance the intensity of the 
beam. The optical thickness normalizes the path length with respect to the 
different scattering and absorption processes that the beam encounters. The 
quantities J and r are functions of both the wavelength of the beam and the 
composition of the medium. 
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Figure 2.3 Graphical illustration of scattering, 
absorption, and emission. 
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The angular distribution of the scattered radiation is described by the phase 
function, P(A;0,£), where 0 is the direction of the incoming radiation, and f is the 
direction of the scattered radiation. It is incorporated into the radiative transfer 
equation through the source function, J. The exact form of the phase function 
depends on both the wavelength and the size and shape of particles of the 
medium. 

The formal solution of the radiative transfer equation is obtained by 
integrating it over the interval (tj, t^), or, 

L(A ;t 2 ;^) = ~ r l> + ~ if hi 

( 2 . 8 ) 

(Chandrasekhar, 1960). The solution consists of two terms. The first term on the 
right hand side represents the direct beam which is attenuated in accordance with 
the optical thickness. The second term represents an augmentation in intensity 
due to external scattering and emission into the beam. The actual solution 
requires knowledge of the behavior of r and J throughout the medium. 

2.2 Radiance Observed bv Nadir-Viewing Satellite 

2.2.1 Scattering and Absorption Mechanisms of & Clear Atmosphere 

It is useful to examine the principal scattering, absorption, and emission 
mechanisms which govern radiative transfer at visible and near-infrared 
wavelengths (0.4 - 1.3 /zm) in the earth's atmosphere under clear-sky conditions. 
Clear atmospheres consist principally of gaseous molecules, particle aerosols (dust), 
water droplets and ice crystals. Gas molecules are the smallest constituent, 
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generally much less than 0.1 fan in diameter. Their spatial and temporal 
distribution depend on the specific gas. Many of the principal gases, such as 
nitrogen, oxygen, and argon, exhibit relatively constant vertical distributions with 
time. Other important gases, such as ozone, carbon dioxide, and sulfur dioxide, 
are variable with time. Water vapor can show significant horizontal and vertical 
spatial variability over time intervals as short as hours (Liou, 1980). 

Aerosols include solid particulates and haze water droplets that range in size 
from 0.1 fan to 1.0 fjcm. The largest atmospheric particles include cloud water 
droplets and ice crystals with mean sizes between 1.0 fan and 10 fan. The 
aerosols are also highly variable in time and space (Iqbal, 1980). 

All the above atmospheric constituents scatter solar radiation. The amount 
and direction of scatter depends on the relative scale of the incident wavelength as 
compared to the size and shape of the particle, as well as the volume density of 
the particles. When the size of the particles is much smaller than the size of the 
wavelength, scattering can be represented by Rayleigh scattering. Most gas 
molecules in the atmosphere contribute to scattering in the visible and 
near-infrared wavelengths approximately according to that mechanism (Slater, 
1980), although for clear skies, the greatest scatter occurs at the shorter 
wavelengths (i.e. 0.4/im). When the particle sizes are comparable to or larger 
than the wavelength, the scattering is termed Mie scattering. Aerosols, cloud 
droplets, and ice crystals scatter solar radiation according to the Mie mechanism 
(Liou, 1980). 

Absorption due to transitions in molecular energy levels occur throughout 
the visible and near-infrared spectrum, although the greatest effect is observed in 
the near-infrared region. Diatomic oxygen has three weak absorption bands 
centered at 0.63, 0.69, and 0.76 fan (Iqbal, 1983). Ozone absorbs from 0.45 to 
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0.77 fan. Water vapor absorbs at 0.72, 0.81, 0.94, 1.10, 1.38, and 1.87 fan. 

Carbon dioxide exhibits weak absorption bands 1.4, 1.6, and 2.0 pan (Liou, 1980). 

Although the integrated effects of all the different types of scattering and 
absorption are a complex phenomenon, the solution to the radiative transfer 
equation is facilitated by treating many of the processes in an additive manner. 

An important quantity to which this applies is the optical thickness. For the 
atmosphere, r(A) can be treated as a sum of optical thicknesses due to Rayleigh 
scattering by gas molecules, r R (A), Mie scattering by aerosols, r M ( A), and 
absorption by gas molecules, r (A) (Liou, 1980), or 

cL 

r(A) = r R (A) + r M ( A) + r ft ( A) (2.9) 

Practically, the total optical thickness of clear skies (low aerosol density and 
water content) is much less than 1.0. For example, observations over southern 
Arizona during the dry season typically range from about 0.5 at 0.4 pan to less 
than 0.1 at 1.5 fan (Slater, 1980), with perturbations occurring at the principal 
absorption bands of water vapor. The decreasing trend is principally due to the 
decrease in Rayleigh scatter with increasing wavelength. Cloudy, aerosol-laden 
skies can possess optical thicknesses greater than 10 (Chahine, 1983). 

2.2.2 Upwelling Radiance at ihg Tod ihfi Atmosphere Under Clear Skies 

The radiance leaving the top of the atmosphere depends on numerous 
complex factors including the solar zenith angle, ground reflectance and 
topography, and atmospheric composition. Despite the complexity, analytical 
solutions for the theoretical nadir radiance observed by a satellite for the case of a 
cloudless sky have been derived by several authors (See for example, Kaufman, 
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1982; Dave, 1980; Otterman and Fraser, 1979, Otterman et al, 1980; and 
Otterman, 1978, 1981). The general approach in most of those models is to 
assume single scattering and to treat the total outgoing radiance as the sum of 
direct and diffuse components. The single scattering assumption has been shown 
to be appropriate for small optical thicknesses (i.e. less than about 0.5) (Bugnolo, 
1960). 

The purpose of this section is not to review the merits of the above 
referenced models, but to select and describe one of them for the purpose of 
providing a framework for incorporation of the stochastic landsurface reflectance 
model presented in Chapter 4. Although any one of the models could be selected, 
the formulation by Otterman and Fraser is presented because it incorporates many 
of the atmospheric processes identified above, as well as the reflectances of the 
target pixel and surrounding pixels. 

Otterman and Fraser (1979) developed an analytically attractive theoretical 
expressions for the nadir radiance as measured by a satellite, for the case when the 
target pixel possesses a different reflectance than the surrounding area. Their 
concern was to account for adjacency effects, that is, the contribution of light 
scattered from areas surrounding of the target pixel into the sensor's field of view. 

The authors' approach was to treat the total nadir radiance observed by the 
satellite as the sum of three components as depicted in Figure 2.4. First, the 
direct beam component, L r (A), represents the portion of solar irradiance reflected 
vertically from the target pixel that is attenuated by atmospheric effects, or 


L r (A) = R(A)E(A) exp[-r(A)]/x (2.10) 
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Target pixel 
reflectance, R(x) 


Figure 2.4 Graphical illustration of direct and scattered diffuse 
components of radiance at the zenith 
(After Otterman and Fraser, 1979). 
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where R(A) is the reflectance of the target pixel, assumed Lambertian, and E(A) is 
the global surface irradiance on the object pixel (Otterman, 1978). 

A second term, LJA), accounts for adjacency effects, or more specifically, 
the portion of diffuse radiance scattered from the surrounding vicinity into the 
column above the object pixel. It is expressed, 


L a ( A ) = [a(A)E(A)/jrr(A)]J ' l-exp[-r(A)/cos4]jcos^[r R (A)P R (A,0 


+ r M^ A ^ P M (A ’^ 2,rsin ^ d £ 


( 2 . 11 ) 


where 

a(A) = average reflectance of the adjacent area, 

P R ,P M = phase functions associated with Rayleigh and aerosol (Mie) 

scattering, respectively, and 
f = zenith reflection angle. 


The size of the area surrounding the target pixel that contributes to this diffuse 
term is discussed by Otterman and Fraser (1979). 

Finally, the third component, L ( j(A), describes the radiance scattered from 
the direct solar beam into the direction of the satellite due to atmospheric effects, 
or 


L d (A) = 


ji{l - exp[- r(A)(l + sec #)]) [ r R U ) P R (A;180' - 0) + r M (A)P M (A;180- - 8)) 
~ r(A)( 1 + /i ) 


( 2 . 12 ) 


where 6 is the solar zenith angle. 
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The total radiance observed by a nadir-viewing satellite is the sum of the 
direct and two diffuse terms, or 

L(A) = L r (A) + L a (A) + L^A) (2.13) 

For the case of optically thin atmospheres, and for clusters of pixels which do not 
exhibit highly contrasting reflectances, the direct beam constitutes the major 
portion of the observed radiance. Exact percentages are difficult to quantify since, 
as indicated by (2.10) through (2.13), their relative proportions are functions of 
numerous variables including surface reflectance, optical thickness, and solar angle. 
Further, most field studies have been conducted over oceans which possess low 
surface reflectance (Chahine, 1983), and thus, they provide limited insight on 
problems over landsurfaces. However, investigations by Otterman and Fraser 
(1979) and Chahine (1983) suggest that the relative magnitude of the direct beam 
ranges from approximately 70 percent of the total nadir radiance for an optical 
thicknesses of about 0.1 and surface reflectance of about 0.4, to as much as 95 
percent for an optical thickness of 0.02 and the same surface reflectance. The 
direct beam percentage decreases with increasing optical thickness and decreasing 
surface reflectance. 


2.3 Characteristics of Landsat Sensors 

Five Landsat satellites were launched between 1972 and 1986 for the purpose 
of investigating earth resources. The technical specifications of the Multispectral 
Scanner System (MSS) and the Thematic Mapper (TM) aboard those satellites are 
well documented (See, for example Freden and Gordon, 1983; Lillesand and Kiefer, 
1987). Several important features relevant to this thesis are summarized below. 


49 



The MSS sensors on Landsats 1, 2, and 3 are line-scanning devices covering 
185 km swaths in four spectral bands. Those include two in the visible spectrum 
at 0.5-0. 6 fan (green) and 0.6-0. 7 fan (red), and two in the near-infrared at 
0.7-0.8 fan and 0. 8-1.1 /zm. The satellite orbits are sun-synchronous, crossing the 
equator at 9:42 a.m. local sun time. The instantaneous field of view (IFOV) is 
square with a ground resolution of 79 meters on a side. The total field of view of 
the scan is 11.6 degrees and the quantization range of the sensor is 64 digital 
numbers (DN). The return period of the sensor for the same area is 18 days. 

The TM sensors on Landsats 4 and 5 have a total of seven bands that 
possess greater sensitivity, finer ground resolution and narrower bandwidths than 
their MSS counterparts. The spectral bandpasses include three in the visible 
region at 0.45-0.52 fan (blue), 0.52-0.60 fan (green), and 0.63-0.69 fi m (red), and 
four in the infrared at 0.76-0.90 /zm, 1.55-1.75 mm, 10.4-12.5 /zm (thermal), and 
2.08-2.35 fan. The TM data are collected at 30 meter ground resolution at a 
quantization of 256 DNs. The orbit is sun-synchronous with an equatorial 
crossing at 11:00 a.m. local sun time, and with a 16 day return period. 
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Chapter 3 

VEGETATION AND SOIL REFLECTANCE 
3.1 Reflectance of Homogeneous Canopies 

The radiation reflected from horizontally homogeneous canopies results from 
the scattering and reflectance properties of the plant components and soil 
background. These properties are both geometric and biophysical in nature and 
thus depend on the species, maturity and overall health of the plant, and on the 
structure and composition of the underlying soil. Geometric plant properties 
include plant architecture, total leaf area, and leaf structure, orientation and 
distribution. Biophysical properties allocate the radiative energy absorbed by the 
plant to important metabolic processes including photosynthesis, respiration, and 
transpiration. Those biophysical properties are manifested in terms of leaf color, 
transparency, temperature, and shape and orientation. Plant properties can vary 
daily and seasonally, in response to soil moisture and nutrients, and to 
meteorologic and climatic conditions. 

3.1.1 Spectral Distribution 

A typical spectral distribution of healthy green vegetation, as compared to 
soil, is shown in Figure 3.1. In the visible spectrum, vegetation reflectance is 
characterized primarily by the absorption of light by chlorophyll in the leaves. 
Strong absorption bands are centered at 0.45 /im (blue) and 0.67 fim (red), 
resulting in a local peak reflectance at about 0.54 fim that corresponds to the 
green portion of the spectrum. In the near-infrared region, vegetation reflectance 
is dominated principally by scattering properties of the internal structure of the 
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Figure 3.1 


1 r- 

MSS red MSS near-infrared 


TM red TM near-infrared 


Typical spectral reflectance functions 
of vegetation and soil. 
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plant, such as the spongy mesophyll of leaves. High scattering in this region 
results in comparatively high overall reflectance by the plant. 

Actual reflectances of plants in the visible and near infrared are well 
documented, especially for crops (Smith, 1983; Myers, 1983). It is typically low in 
the visible spectrum (< 30%) and higher in the near infrared region (40 - 50%). 
Canopy reflectances are generally much lower than those measured for individual 
leaves (Dickinson, 1983). Kondratyev (1969), Gates (1980), Ross (1981) and Iqbal 
(1983) provide summaries of reflectances for various vegetation types. 

In contrast to vegetation, soil reflectance generally exhibits increasing 
reflectance with wavelength, as shown in Figure 3.1. Overall magnitudes are 
governed by grain size distribution, soil moisture, organic content, etc. (Myers, 
1983). Some of the important factors are examined later in this chapter. 

3.1.2 Reflectance Models 

Numerous radiative transfer models for horizontally homogeneous canopies 
have been developed in terms of various plant properties and background soil 
reflectance. Typically, homogeneous canopies have been modeled as a diffusing 
medium with absorbing and scattering properties. Excellent reviews of these 
models are provided by Smith (1983) and Ross (1981). Suits (1972) envisioned a 
plant canopy as an infinitely extended plane-parallel medium with homogeneous 
geometric properties. Verhoeff and Bunnik (1981) extended the Suits model to 
include the effect of leaf angle distribution. Dickinson (1983) applied the 
two-stream approximation for radiation transfer in the atmosphere (Meador and 
Weaver, 1980) to plant canopies employing the leaf area index as a measure of the 
optical depth. Recent literature has increased the sophistication of those earlier 
models to include the modeling of bidirectional reflectance (Camillo, 1987; 
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Walthall et al., 1985; Chen, 1985; Simmer and Gerstl, 1985; Vanderbilt and 
Grant, 1985, Strebel et al., 1985; Gerstl and Zardecki, 1985; and Reyna and 
Badhwar, 1985). 

Semiempirical formulas for the total radiation fluxes have been proposed for 
practical applications. The attenuation of radiation as it passes through a plant 
stand has been typically described in terms of some form of Bouguer's Law such 
as that proposed by Monsi and Saeki (1953), or 

E(A;C) = E q (A) exp[-p(A)C] (3.1) 


where 

E 0 (A) = intensity of incoming radiation at top of canopy 

E (A;C) — intensity of radiation at a penetration level associated with leaf 
area index, ( 

Ai(A) = experimentally determined extinction coefficient 

When the upwelling radiance from the canopy represents only the backscattered 
solar radiation from the plant (i.e. no soil effects), then the plant reflectance, 
R m (A;C)i is defined 


R m (A;0 =l-e^< A >< (3.2) 

Other semiempirical formulas which account for the partitioning of transmitted, 
scattered, and absorbed radiation, have been proposed (Ross, 1981). 

Attention has also been focused on the invertibility of reflectance models of 
homogeneous canopies for estimating plant parameters such as leaf area index, 
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biomass, and leaf angle distribution (Goel and Strebel, 1983; Goel et al., 1984; 

Goel and Thompson, 1984a, 1984b, 1984c; Lang et al., 1985). Goel and Thompson 
(1984c) have shown that such parameters can, in principle, be estimated using 
multispectral data at approximately 25 combinations of solar zenith angle and 
viewing angle. 


3.2 Reflectance of Discontinuous Canopies 

Most natural landscapes will vary both horizontally and vertically in species 
and/or vegetation density. Modeling this situation has received considerably less 
attention than that for homogeneous canopies. Statistical techniques have been 
employed for classifying landscapes based on their similar spectral signatures (See 
summary by Lillesand and Kiefer, 1987). However, such methods require the 
identification of training samples and therefore can not be adopted for natural 
landscapes in which all target pixels possess unique spectral characteristics. 

Radiative transfer models for non-homogeneous canopies have been developed 
by extending homogeneous canopy models and including three-dimensional 
scattering functions (Ross, 1981; Kimes et al, 1985). Such models can be solved in 
a few cases where plant distribution is periodic, such as for sown crops (Suits, 
1985). Inversion of such models, using only multispectral observations, also 
require a large number of data at different angles (Goel and Thompson, 1984c). 

For regional scale investigations of natural systems, the application of the 
above radiative transfer models, using Landsat data or any other current satellite 
multispectral sensor (i.e SPOT, AVHRR), is clearly infeasible. The first 
constraint is simply that any of those systems provides only one observation per 
overpass. Additionally, however, the regional-scale heterogeneity in plant species, 
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size and architecture would prohibit making any assumptions on the distribution 
of plant parameters. 

3.2.1 Proportion Models 

In many regional scale hydrologic investigations, it is generally sufficient to 
know only the amount of vegetation and the outgoing radiance at its surface, and 
not the complex architecture within the plant. In such cases the pixel has been 
described in terms of its bulk components, that is vegetation, soil, and shadow. 

A simple approach has been to assume that the surface reflectances of 
individual plant clusters and soils within a pixel are constant, and that the total 
spectral response is a linear combination of the individual spectral responses of its 
components (Horwitz et al., 1971; Nalepka et al., 1972; Work, 1974; McCloy, 1980; 
Dozier, 1981; Ungar and Bryant, 1981; Chhikara, 1984, Lenningtion et al, 1984). 
The total radiance emitted from a pixel, M(A), containing n cover types can be 
expressed 


n 

M(A) = l fjMjfA) (3.3) 

i=l 

where 

fj = fraction of area covered by type i, and 
Mj(A) = radiance emitted from cover type i in band A. 

For the simple case of vegetation and soil cover, Equation (3.3) becomes: 


M(A) = mM m (A) + (1 - m)Mg(A) (3.4) 
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ni = percent vegetation cover 

M m (A) ,M (A) = radiances emitted from vegetation and soil, 

respectively, in band A. 

Perhaps the earliest development of proportion estimation can be attributed 
to Horwitz et al. (1971) and Nalepka et al. (1977) who also termed the method 
"mixtures estimation". One of its first applications was in identifying subpixel 
scale ponding and wet marshes in glaciated prairies (Work, 1974). 

McCloy (1980) later proposed that under conditions of negligible canopy 
transmission or multiple reflection, the response proportions of the various land 
covers will closely approximate the physical proportions of each type of cover. He 
suggested that up to four sub-pixel categories be used including three levels of 
vegetation greenness cover and one soil background cover. Ungar et al. (1981) 
reported limited success delineating forest canopy types in Maine using a similar 
approach which they termed the "Fanning algorithm". The fractional area was 
determined by minimizing the error between observed and theoretical radiances. 

3.2.2 Geometric Models 

In an extension to the linear proportion model described above, some 
investigators have considered the shadow cast by vegetation as an additional 
component to the total radiance. These models abstract clumps of vegetation as 
three-dimensional geometric shapes on flat horizontal surfaces with constant 
reflectances (Otterman, 1981, 1984; Otterman and Weiss, 1984; Strahler and Li, 
1981; Li and Strahler, 1985). The distribution of the elements themselves can be 
geometric, as in the case of row crops, or statistical, as for natural vegetated 
landscapes. In most cases, the surfaces are assumed Lambertian and scattering 
between vegetation and soil is neglected. 
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Geometric models have been successfully used to describe much of the 
variability of semivegetated landscapes by altering the shape and density of the 
geometric elements. Otterman (1981, 1984) and Otterman and Weiss (1984) 
envisioned forests and desert vegetation as thin vertical cylinders. They accounted 
for the shadowing effects of vegetation, but assumed that the horizontal area of 
the plants was negligible. That model is thus not directly applicable to the 
determination of fractional vegetation cover. 

Richardson et al (1975) modeled crop canopies as rectangular rows, 
neglecting scattering between the crop and soil. Jackson et al (1979) extended the 
above model to include shadowed canopy effects. 

Strahler and Li (1981) and Li and Strahler (1985) modeled conifer forests as 
randomly located cones of similar shape and random height. They determined 
from simple geometry the shadow cast by the cones on the soil background or on 
other cones. The total radiance emitted by a pixel was assumed to consist of four 
components: illuminated background, illuminated cones, shadowed background and 
shadowed cones. Vegetation parameters including percent cover and average tree 
height were then estimated using assumed values of component reflectances. 

3.3 Spatial Pattern of Vegetation 

The application of geometric models to natural watersheds generally requires 
assumptions on the statistical distribution of plant spacing. Several authors (see, 
for example, Whittaker, 1975; Diggle, 1983) have focused on the problem of fitting 
stochastic models to spatial point patterns of natural (undisturbed) vegetation. 
Diggle (1983) found that the most appropriate stochastic representation for a given 
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situation was species dependent, with several forest species following a Poisson 
distribution or a Poisson cluster process. 

Significant progress in incorporating statistical spatial distributions into the 
analysis of remotely sensed data has only recently been achieved (Strahler and Li, 
1981; Li and Strahler, 1985; Woodcock, 1985). Li and Strahler (1985) and 
Strahler and Li (1979) assumed a homogeneous Poisson distribution of conifer tree 
locations. Woodcock (1985) used a similar model to examine the relation among 
the scale of pixel components, resolution size, and two indicators of spatial 
correlation: the variogram and the local vfl.ria.nrp 

3.4 S2ii Reflectance 

Soil background reflectance often exhibits high variability in semivegetated 
scenes. It varies over a wide range of length scales due to changes in the physical 
structure and chemical composition of the near surface soil. Small-scale 
perturbations (meters or less) occur with changes in mineral content, water and 
organic matter content, particle size, soil texture and surface roughness. 

Numerous experimental investigations have examined the relationship between bare 
soil reflectance and those parameters (for a summary, see Myers, 1983). 

Soil reflectance also varies at larger geophysical scales. For example, the 
presence of hills will change soil reflectance due to an effective altering of the 
illumination and viewing angles. Changes in elevation, slope, and aspect will 
influence soil moisture and organic matter content. Subsurface variations in 
geologic formations will affect mineral content. Those geophysical factors impose 
a spatial correlation to soil reflectance at scales of 10 to 10^ meters. 
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In addition to spatial correlations, soil reflectance can also be 
cross - correlated at different wavelengths. These wavelength dependent correlations 
are manifested in multispectral scattergrams of real scenes through the preferred 
location and orientation of bare soil pixels (see, for example, Kauth and Thomas, 
1976; Richardson and Wiegand, 1977; Huete et al, 1985). For red-infrared 
scattergrams of typical semivegetated scenes, the data often take on the form of a 
triangle whose base consists of a straight line emanating from approximately the 
origin. Researchers have identified that base line, consisting primarily of pixels 
containing bare soil and dry vegetation, as the "soil line". 

Analysis of previously published theoretical and experimental studies (for 
example Cierniewski, 1987; Smith, 1983; Bowers and Hanks, 1965; Skidmore et al, 
1975) indicates that for a given type of soil variability, the soil reflectance at one 
wavelength is often functionally related to the reflectance in another wavelength. 

In many cases, the relation can be approximated by a simple linear expression 

R(A 2 ) = aR(Aj) + 7 (3.5) 

where the slope, o, and intercept, 7 , are coefficients dependent on both the wave- 
length and the type of variability. Thus, variability of any one soil parameter can 
lead to a representative "line" in a two-dimensional scattergram. 

For instance, Figure 3.2 contains three hypothetical visible— infrared scatter- 
grams, representing three different scenes, in each of which only one soil parameter 
is allowed to vary. In Scene 1 only the amounts of two minerals are allowed to 
vary, while all other soil parameters such as surface roughness, moisture, etc. are 
held constant. The resulting scattergram forms a "soil mineral line" in which the 
end points approach the respective reflectances of the pure minerals. The shape 
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and orientation of the line may be linear (as drawn) or nonlinear, and are 
determined by the location of the end points and the nature of the mixing of the 
two mineral types (Smith et al, 1985). Pixels lying between the end points will 
contain mineral amounts proportional to their distance along the line. 

Scene 2 contains hypothetical pixels in which only soil moisture is allowed to 
vary. Soil moisture increases the radiation absorption capacity of the soil in the 
visible and near infrared regions. Analysis of published experimental data (Bowers 
and Hanks, 1965; Skidmore et al, 1975) indicates that, for many soil types, 
equation (3.5) is applicable. Thus the locus of points for Scene 2 pixels should 
form a "soil moisture line" as indicated in Figure 3.2 with the pixels along the left 
portion of the line containing higher soil moisture than those to the right. 

Finally, Scene 3 contains pixels in which only surface roughness varies. Soil 
reflectance generally decreases with increased surface roughness due to the increase 
in shadow cast by the soil particles and aggregates onto itself (Myers, 1983, 
Cierniewski, 1987). The resulting "soil shadow line" is approximately linear with 
an intercept of near zero. The linearity occurs since the amount of shadow caused 
by the soil aggregates is practically the same for the range of wavelengths being 
considered. In fact, a near zero intercept for straight soil lines of real scenes (with 
low diffuse radiation) may be an indication that soil shadow induced by its 
physical structure is the dominant source of soil reflectance variability. 

In actuality, real soil scenes contain a composite of several types of varia- 
bility. The corresponding soil line is generally linear in the mean although 
considerable scatter can exist (Kauth and Thomas, 1976; Richardson and Wiegand, 
1977) A unique soil line will exist only if either 1) one dominant type of soil 
variability is occurring or 2) the scatter due to the different types of soil 
variability act in the same direction. 
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3.5 Empirical Indicators of Vegetation Amount 


Numerous vegetation indices have been proposed in recent years as 
qualitative indicators of green vegetation. The purpose has been to reduce the 
several multispectral bands to one value to estimate vegetation parameters such as 
biomass, leaf area index, or percent cover. The relationship among vegetation 
indices and the plant physical properties was investigated by Choudhury (1987) 
using a two-stream approximation. The effect of soil background reflectance was 
examined by Huete (1988). Perry and Lautenschlager (1984) summarize the many 
different vegetation indices and describe their relationship to each other. Three 
such indices are the normalized vegetation index, the perpendicular vegetation 
index, and the Kauth-Thomas Greenness index. 

3.5.1 Normalized Difference Vegetation Index 

Of the many indices proposed, the normalized vegetation index (NDVI) has 
evolved as a practical popular choice for use in regression with vegetation 
parameters. It is of the form 

R(NIR) - R(VIS) 

NDVI = R(NIR) + R( VIS ) 

where R(NIR) and R(VIS) represent the pixel reflectances or the digital numbers 
(DNs) in the near-infrared and visible (red) ranges, respectively. Low NDVI 
indicates low vegetation amount, whereas high NDVI indicates either high 
vegetation amount or high productivity (Curran, 1980). Sellers (1985) discussed 
the functional relationship between the NDVI and several vegetation parameters, 
including the leaf area index. Tucker et al. (1983) correlated the NDVI 
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(computed from NOAA's AVHRR data) to actual biomass obtained from field 
sampling in a semi-arid region of Senegal, West Africa. The effect of soil 
background on NDVIs computed from hand-held or airborne radiometer data was 
investigated by Elvidge and Lyon (1985) and Huete et al. (1985). 

3.5.2 Perpendicular Vegetation Index 

Richardson and Wiegand (1977) proposed the perpendicular vegetation index 
(PVI) as a measure of plant development. Application of this index first requires 
the establishment of a background soil line by linear regression of MSS bands 2 
and 4 using bare soil pixels. The soil line is thus a straight line stemming from 
near the origin. The PVI is the perpendicular distance from the soil line to the 
actual data point which contains vegetation, and is defined, 


PVI = [( Rg2 - R p2 ) 2 + (R ?4 - Rp 4 )Y /2 (3.7) 

where ^g2’^g4 = re ^ ectances of soil background in bands 2 and 4, 

respectively, corresponding to the data point. 

R p2 ,R p4 = reflectances of data point in bands 2 and 4, 

respectively, perpendicular to R _ and R„, on the soil 

g 2 g 4 

line. 

Richardson and Wiegand (1977) regressed PVI with percent cover of 
sorghum with a correlation coefficient of 0.57. Theis et al. (1984) studied the 
effect of vegetation and soil moisture on PVI. Rosenthal et al. (1985) recently 
used the PVI to investigate crop height and biomass. 
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3.5.3 Kauth-Thomas Greenness Index 

Kauth and Thomas (1976) applied Gram-Schmidt orthogonalization to the 
original four Landsat bands resulting in a new four-dimensional space termed 
"tasseled cap". The procedure, which is similar to principal components except in 
the order of calculations, essentially rotates the data so that most of the 
variability can be explained in terms of four indices: greenness (GI), brightness 
(BI), yellowness (YI), and nonsuch (NI). The first two of these indices are 
defined, 


BI = 0.332 DN1 + 0.603 DN2 + 0.676 DN3 + 0.263 DN4 (3.8) 

GI =-0.283 DN1 - 0.660 DN2 + 0.577 DN3 + 0.388 DN4 (3.9) 

where DN1, DN2, DN3, and DN4 represent the digital counts of the three visible 
and one near infrared bands of the MSS scanner, respectively. A similar set of 
equations has been developed for the Thematic Mapper (Crist, 1983; Crist and 
Cicone, 1984). 

The Kauth-Thomas Transformation has been used by numerous investigators 
to model various crop parameters including crop development, moisture stress, 
yield and crop classification (Ezra et al., 1984). Huete et al. (1984, 1985) in a 
series of small scale experiments of wooden boxes filled with soil, showed high 
correlation of GI with percent cover. Musick (1984), however, using Landsat MSS 
data over New Mexico, was unable to achieve consistent differentiation between 
arid rangeland cover changes using GI. 
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Chapter 4 


A STOCHASTIC-GEOMETRIC LANDSURFACE REFLECTANCE MODEL 
FOR SATELLITE REMOTE SENSING 


This chapter describes a stochastic-geometric landsurface reflectance model 
that can be incorporated into the theoretical nadir radiance model described in 
Section 2.2.2. The goal is to provide a framework for investigating the conditions 
under which subpixel variability is represented in satellite-observed radiance. The 
first part of this chapter presents a stochastic canopy-soil reflectance model for 
describing regional landsurface variability. The reflectance model is then coupled 
to the radiance model through the reflectance term. 



Many semivegetated landscapes are characterized not only by variations in 
plant size and density, but also by complex variations in the physical properties 
controlling the reflectance of the plant-soil medium. The complexity arises 
principally from the highly random nature of many subpixel scale properties of the 
plant (e.g. species, plant geometry, leaf area, shape, and orientation) and the soil 
(e.g. surface roughness, texture, organic matter content, and moisture content). 
Such variations, induced by changes in topography and climate, can have a 
significant influence on the interpretation of scenes and therefore must be 
recognized when applying canopy reflectance models to regional scale problems. 

One approach for accommodating the random subpixel fluctuations in plant 
and soil properties, while keeping the number of model parameters to a minimum, 
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is to conceptualize the semivegetated landscape as a stochastic-geometric reflecting 
surface. This approach is similar to the geometric models described in Section 

3.2.1 in that the landscape is represented by the bulk plant and soil properties 
that dominate scene reflectance, that is, fractional cover, plant geometry, and 
plant and soil reflectance. However, it extends those models by treating any one 
of the bulk properties as a potential random variate. Treating the properties as 
random variables provides a flexible means of characterizing the scene without 
having to prescribe an inordinate number of detailed vegetation and soil 
parameters (e.g., leaf area and orientation, soil roughness, etc.). Such an approach 
is practical in many mesoscale hydrologic investigations in which detailed 
description of the surface is not necessary, nor is it feasible with current satellite 
sensors. 

4.1.1 Plant and Shadow Geometry 

The geometric aspects of the reflectance model include both the plant 
geometry and the spatial distribution of the plants. Since the plants are 
represented as solid three-dimensional objects, they can cast shadows onto 
themselves, onto the soil, and onto adjacent plants, according to their shape and 
spatial distribution, and to the solar angle. 

The four principal reflecting surfaces of the model are the illuminated green 
canopy, the shadowed green canopy, the illuminated soil background, and the 
shadowed background. All surfaces are assumed to be Lambertian and scattering 
among geometric elements is neglected. The term "soil background" denotes that 
the landsurface between the green plants includes not only bare soil, but also a 
mixture of sparse grasses and senesced vegetation. 
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The geometry of the plant canopy and its shadow is shown in Figure 4.1 for 

the particular case of cones. Li and Strahler (1985) has shown that, assuming a 

plant diameter D, and height H, the portions of the illuminated plant, A. , and 

T 

shadowed plant, , are respectively, 


A. = [*/2-*]D 2 /4 (4.1) 

A, = [*/2 + x)D 2 /4 (4.2) 

s 

where 

\ = sin _1 [D/Htan0| (4.3) 

Geometric relationships for other shapes, such as cylinders and spheres, are 
provided in Section 6.2.3 and Table 6.1. 

4.1.2 Reflectance of a Pixel 

Depending on its size, an individual pixel can possess any number of plants 
and shadows, or fractions of plants and shadows. The total reflectance of an 
individual pixel, R(A,x), is assumed to consist of an area weighted linear 
combination of the average reflectances of the four landscape components. 
Mathematically, 

R(A,x) = E fj(x)Rj(A,x) (4.4) 

where, 

A = wavelength, 

x = spatial coordinate of the center of the pixel, 

Rj(A,x) = average bulk reflectance of cover type i in pixel centered at x. 
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Illuminated soil background 

Shadowed soil background 


Figure 4.1 Geometry of canopy-soil reflectance model 
for the case of cones. 
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fj(x) — fraction of pixel centered at x possessing the cover type i as 
follows 

i = 1 designates illuminated green vegetation, 

i = 2 designates shadowed green vegetation, 

i = 3 designates illuminated soil background, and, 

i = 4 designates soil background shadowed by green 

vegetation. 

The percent covers are constrained by, 

£ f jU) = 1 (4.5) 

Average bulk reflectances and percent covers are defined as follows, 



where, 

r j( A,ii) = reflectance of point u, given it possesses cover type i, 

Aj(x) = total area of the pixel with cover type i, and 
Ap = total area of a pixel, 

Equations (4.4) through (4.7) are applicable for any pixel in the scene. 
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4.1.3 Smtial and Spectral Distribution of Scene Variables 

The equations representing total pixel reflectance, (4.4) through (4.7), are 
completely general in that they do not specify whether the scene variables are 
deterministic or random. The quantities f.(*) and Rj(A,x) can thus represent 
constants, deterministic variables, or random variables depending on the nature of 
the scene. For instance, the shape and spatial distribution of agricultural plants 
can often be prescribed by a regular geometry, such as rectangular rows. In 
natural areas, however, plant height and regional distribution can more 
appropriately be represented by random functions as noted in Section 3.3. Plant 
and soil reflectance may also be deterministic or random as noted above. 

When one or more of the scene variables are random, the reflectance model 
takes on a stochastic nature in that the values of fj(x) and Rj(A,x) can become 
correlated in a number of possible ways. For instance, there may exist i) spatial 
correlations in plant and soil reflectance, ii) cross-spectral correlations between 
plant reflectance and soil reflectance, and iii) cross-correlations between fractional 
covers (i.e. vegetation and shadow). Some of those correlations are examined 
further in Chapters 6 through 8. The stochastic nature of the reflectance model 
thus provides a mechanism for solving the inverse problem as described later in 
Chapters 6 and 7. 

4.2 Coupled Landsurface-Atmosphere Radiative Transfer Model 

The landsurface reflectance model described in (4.4) is coupled to the 
atmospheric radiance model through the target reflectance parameter, R(A), in the 
direct beam equation (2.10), and through the adjacent area reflectance parameter, 
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a(A), in the diffuse scattered equation (2.11). For instance, inserting (4.4) into 
(2.10) yields, 

L.(A,x) = E fj(x)Rj(A,x)E(A)exp[-7-(A )]/ r (4.8) 

The diffuse radiance scattered from the surrounding vicinity can be written, 

L a (A,x) = Efj(x)a i (A,x)E(A)/xr(A) J ^ Jl-exp[-r(A)/cos £]Jcos £k R P R (A,f) + 

1 0 

r M( A ) P M (A ’^ 2TSin ^ d £ (4.9) 

where fj(x) represents the average amount of cover type i in the area surrounding 
the target pixel located at x, and a^(A,x) represents the average reflectance of 
cover type i in that same area. The backscattered solar radiance, L d (A), described 
in (2.12), contains no landsurface parameters and is thus unaffected by the nature 
of the landscape. 

Since the focus of this report is on the landsurface cover, it is useful to 
separate the landsurface properties from those of the atmosphere. Equations (4.8) 
and (4.9) can be rewritten, 

L r (A,x) = (Sf.(i)R i (A lJ £)]L;[A,E,r] (4.10) 

L & ( A ix) = S fj(x)aj(A,x) L a [A,E,r,P^,P j^,£] (4-11) 

where LJ. and represent the expressions in (2.10) and (2.11), respectively, 
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without the reflectance terms. The total radiance observed by the satellite is 
thus, 



4.3 Eff§£t Qf Subpixel Variability Qn Observed Radiance 
Equations (4.10) through (4.12) can be used, at least in a qualitative 
manner, to examine the theoretical relationship between fractional cover and 
satellite-observed radiance, given the limiting assumptions of both the landscape 
reflectance model and atmospheric radiance model. For instance, (4.10) indicates 
that the direct beam, the dominant component of the observed radiance for 
optically thin atmospheres, is linearly proportional to the amount of the individual 
fractional covers of the target pixel (i.e. the fAs). However, that proportionality 
does not exist when the adjacency effects of (4.11) are included. The magnitude 
of adjacency effects will depend on both the amount of atmospheric scattering and 
on the relative contrast between the target pixel reflectance and the surrounding 
area reflectance (Otterman and Fraser, 1979). 

Since the diffuse terms are small compared to the direct beam for clear skies, 
it can be argued from (4.12) that adjacency effects will not be important when the 
average reflectance of the surrounding pixels approximately equals the total 
reflectance of the target pixel, even though the distribution of ground cover in any 
cluster of pixels may differ. In such cases, 

R(A,x) ~ a(A,x) (4.13) 

and the total observed radiance can be expressed, 
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Since the satellite data provides observations in terms of DN values, equation 
(4.14) can be rewritten in the form, 


DN(A,x) 



(4.15) 


where 1^.(A) and 1^(A) are calibration coefficients for a particular horizontally 
homogeneous atmosphere and satellite sensor. Further, since one is primarily 
interested in the landscape fractional cover, then (4.15) can be expressed, 


DN(A,x) = E fj(x)Rj(A,x) + 1^(A) (4.16) 

where Rj(A) represents an "effective" reflectance, that is, the reflectance of 
component i multiplied by the coefficient 1^.(A). Equation (4.16) implies that for 
horizontally-homogeneous, optically-thin atmospheres, the satellite observed 
radiance is linearly proportional to subpixel fractional cover when the landscape 
reflectance does not possess sharp contrasts. That preliminary conclusion has 
important implications for the interpretation of scattergrams obtained from 
satellite data as discussed later in Chapter 8. 




Chapter 5 

SIMULATION OF RED-INFRARED SCATTERGRAMS 
OF SEMIVEGETATED LANDSCAPES 


One application of the stochastic-geometric canopy reflectance model is the 
investigation of the structure, or physical basis, of red-infrared scattergrams of 
semivegetated landscapes. That is achieved by using the reflectance model 
presented in Section 4.1 with typical values to simulate different scenes, and then 
comparing the shape and common features of the corresponding red— infrared 
scattergrams. An understanding of the influence of a given random variable is 

obtained by altering one of its statistics (e.g., variance), while holding all others 
constant. 

The following section presents the results of five different simulations in 
terms of the nadir reflectance model. Atmospheric effects are not considered and 
it is assumed that all areas of the scene are equally illuminated. The input values 
of the different reflectance variables are provided in Table 5.1. While those scenes 
represent only a few selected scenarios, they were chosen to demonstrate the 
important general relationship between the major landscape variables and their 
effect on the scattergram. 

Scenes are generated as follows. A scene consists of eight segments, each 150 
meters square with one meter square grids. Each segment within a scene is 
assigned an identical soil background reflectance distribution. Next, trees of fixed 
height (3.5 m), represented by square cylinders, are superposed on the soil 
background of each segment according to a Poisson distribution having a different 
arrival rate for each segment. The area of the canopy (or cylinder) is fixed at 
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1 m . The shadowed portions of the scene are then determined based on the 
geometry of the tree and the solar zenith angle (9 — 30 degrees). Typical values 
of soil and vegetation reflectance are assumed (Ross, 1981; Myers, 1983; Smith, 
1983) as indicated in Table 5.1. All surfaces are assumed to exhibit Lambertian 
properties. A graphical sketch of a typical pixel in the simulated scene is shown 
in Figure 5.1. 

Once the scene is generated at the 1 meter grid level, it is then aggregated 
to pixel sizes of 5, 10, and 30 meters square by averaging the reflectance values of 
the components of the grid. A comparison of the pixel scales for the different 
levels of aggregation is shown in Figure 5.2. The latter two scales represent 
SPOT and Thematic Mapper satellite pixels, respectively. Information regarding 
subpixel variables is recorded at each level of aggregation. The computer code for 
the generation of scenes is provided in Appendix D. 

The simulations are based on the reflectance equations described in (4.4) 

through (4.7). To facilitate the explanation of the individual cases, the notation is 
changed as follows, 

nij = fraction of illuminated vegetation cover in pixel centered at x, 

m s = fr ac ti° n of shadowed vegetation cover in pixel at x, 

Sg = fraction of soil shadowed by vegetation cover in pixel at x, 

gj = fraction of illuminated soil background in pixel at x, 

V(A*) = average bulk reflectance of illuminated and 

shadowed vegetation, respectively, in band A, in pixel 
centered at x, 

R e (A,x), R (A,x) = average bulk reflectance of illuminated and 

6 i ®s 

shadowed soil background, respectively, in band A, in 
pixel centered at x. 


77 




78 




E 

o 

CO 


c 

o 

cn 

O' 

o 

k- 

O' 

O' 

03 

O 

CO 

"5 

> 

<D 



79 


Figure 5.2 Comparison of pixel scales for different 
aggregations of simulated scenes. 



The reflectance equation of an individual pixel given in (4.4) can be written 
in terms of the new notation as, 



The average bulk reflectances axe computed, for example, in the case of soil 
reflectance, 


R gj( A ’-) = n I r gl ( A ’U)%) (5.2) 

n 


where, 

n 


r g ( A >u) 
g I 

%(u) 


= the number of grids in the pixel centered at x occupied by 
illuminated soil, 

= illuminated soil reflectance at point n, and 

= delta function equal to 1, if point u is occupied by illuminated 
soil, or 0 if it is not. 


When the bulk reflectances are not spatially variable, that is, when they are only 
a function of wavelength, then for clarity, they will be written without the 
coordinate variable x. For example, R (A, 2 c) becomes R (A) when soil reflectance 

6j gj 

is not spatially variable. 

The fractional cover of an individual pixel is defined, 

A 8, (X) 

*1 = < 5 ' 3 > 

P 
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This case represents an idealized two-component situation in which the 
vegetation and soil each have a constant reflectance over the entire scene, and 
observations are from the nadir. The sun is near zenith resulting in no shadows in 
the field of view. Hence, the only random variable is percent cover. The equation 
expressing total reflectance from an individual pixel is taken from equations (5.1) 
through (5.7) with m = g = 0, or 

b S 



The red-infrared scattergrams for the Case I simulation are shown in Figures 
5.3-a,b,c for the levels 5, 10, and 30 aggregation. They indicate that all the data 
points fall on a straight line whose end points represent pixels containing the 
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Figure 5.3-a,b,c Red-infrared scattergrams, Case I simulation: variable 

percent cover, constant soil and vegetation reflectances; no 
shadows. 



maximum (upper left) and minimum (lower right) percent vegetation cover within 
the scene. The length of the line decreases with increasing aggregation since the 
standard deviation of the canopy density decreases as the pixel size increases. The 
percent cover of any pixel lying along the line is proportional to the distance 
between the end points. 

5.2 Case II - Variable Soil Reflectance. No Shadows 

In addition to changing vegetation cover, the Case II simulation includes the 
effect of spatial variability of soil reflectance. Both small scale (subpixel) and 
large scale (regional) variations are incorporated by treating soil reflectance as a 
two-dimensional random field with a prescribed covariance structure. 

While various functional forms might be applicable, for demonstration 
purposes, the Case II simulation assumes that soil reflectance is normally 
distributed with an exponential covariance function. It is expressed 

cov(d) = a 1 exp(-/?d) (5.9) 

where a 2 = the variance of the soil reflectance distribution, 

(3 = inverse length scale of the covariance function, and 

d = distance between two points in the scene. 

The simulated bare soil scene for the red band is shown in Figure 5.4. That 
scene, generated using the Turning Bands Method (Mantaglou and Wilson, 1982), 
contains a mean reflectance (0.15), standard deviation (0.023), and exponential 
covariance. A similar scene (not shown) was generated for the infrared band. 
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Figure 5.4 Hypothetical segment of bare soil scene, red band. 
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The mean and standard deviations of both scenes are based on the hypothetical 
soil reflectance curve shown in Figure 5.5, which indicates a typical range of 
reflectances for a soil with variable properties (e.g., soil moisture or surface 
roughness) (Myers, 1983). 

The length scale of the exponential covariance function was chosen to be 20 
meters. While that might represent some geophysical scale, for the present case, 
it was chosen for convenience. It is much larger than the grid scale of one meter, 
and the two smaller aggregations (5 and 10 meters), but smaller than the largest 
aggregation (30 meters). Thus, the choice of that scale allows one to examine the 
relation between covariance length scale and pixel size. 

The total reflectance from a given pixel in the Case II simulation is, 



The results of the Case II simulation for all three aggregations are shown in 
Figures 5.6-a,b,c. (Regular spaces in the scattergrams, especially at lower level 
aggregations, are due to finite increments in percent cover as limited by the level 
of aggregation. This effect occurs in subsequent cases as well.) They indicate 
that the red-infrared scattergram tends to take on the characteristic shape of a 
triangle. The top of the triangle represents full canopy cover, and the base 
represents the minimum vegetation cover in the scene. For areas in which it can 
be assumed that bare soil exists, the base of the triangle represents the classic 
"soil line." 

While Case II is still a relatively simple example, it demonstrates the 
usefulness of the scattergram for explaining subpixel variability. For instance, all 
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Figure 5.6--a.b,c Red -infrared scattergrams, Case II simulation: variable 

percent cover and soil reflectance, constant vegetation 
reflectance; no shadows. 



pixels falling on a line parallel to the soil line will have the same percent 
vegetation. A second observation is that all pixels falling on a straight line 
extending from the top of the triangle to the soil line will have the same value of 
average soil reflectance. The above interpretations of the scattergram are 
indicated on Figure 5.7 (an expanded version of 5.6b) for the level 10 aggregation. 

The importance of pixel scale relative to the covariance function is seen in 
the size of the triangles at different levels of aggregation. The scattergrams 
indicate that the length of the soil line and, hence, the width of the triangle 
decrease with increasing aggregation. Both the standard deviation and the 
covariance length scale contribute to that effect. Since scenes composed of large 
pixels average over a greater area than scenes with small pixels, statistically, one 
can expect the former case to have a lower standard deviation. However, that 
effect is mitigated by the covariance length scale. Scenes with small length scales 
(relative to pixel size) will exhibit short soil lines, while scenes with large length 
scales will exhibit long soil lines. 



In addition to variable percent cover and variable soil reflectance, the Case 
III simulation introduces variable vegetation reflectance and examines its effect on 
the red -infrared scattergram. Vegetation reflectance will change at small and 
large spatial scales due to variations in a number of plant parameters, including 
plant species, leaf reflectance, growth stage, plant architecture, leaf orientation and 
distribution, leaf area, and plant stress (Ross, 1981). Regional scale variations in 
the pattern of natural vegetation and dominant species are influenced by elevation, 
gradient, and local climate (Whittaker, 1975). 
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RED REFLECTANCE (%) 


Figure 5.7 Interpretation of scattergram, Case II simulation 
level 10 aggregation: variable percent cover and 
soil reflectance, constant vegetation reflectance: 
no shadow. 
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As for soil, Case III treats the variation in vegetation reflectance as a 
normally distributed random variable with an exponential covariance structure. It 
further assumes that reflectances in the infrared and red bands are linearly related 
with negative slope. That relationship is not intended to represent all types of 
vegetation variability, but may be a simple approximation for some cases. For 
instance, increases in leaf area are generally associated with decreases in red 
reflectance and increases in infrared reflectance (see for example, Colwell, 1974; 
Hall, 1984). 

For Case III, the total reflectance of a given pixel becomes 

R(A,x) = mR (A,x) + (1 - m)R (A,x) (5.11) 

*i h 

where the three random variables are percent cover (m), vegetation reflectance 
(R m ) and soil reflectance (R^ ). 

The scattergram for Case III is presented in Figures 5.8— a,b,c for all levels of 
aggregation. The difference from Case II is that the top of the triangle has spread 
open, resulting in a quadrilateral data plot. An envelope curve along the top of 
the quadrilateral represents pixels of maximum vegetation cover. For scenes 
containing full canopy cover, that locus of points can be considered the "canopy 
line" analogous to the soil line at the base of the quadrilateral. 

It is noted that for all three non shadowed cases (I, II, and III), neither 
plant geometry nor spatial distribution play a role in the shape of the scattergram 
or the relative location of a given pixel. Similar scattergrams could have been 
achieved using any plant geometry (e.g., spheres or cones) or spatial distribution 
(e.g., row crops with any orientation) as long as the distribution of the 
reflectances and the percentage of vegetation cover were the same. 
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5.4 Case IV - Shadowed Soil Background. Constant Vegetation and Soil 

Reflectance 

Shadows cast by vegetation can be an important component of total pixel 
reflectance. Shadows change diurnally with the position of the sun and with the 
amount of diffuse solar radiation. Important seasonal changes occur both with the 
sun's migration and with changes in plant structure. 

Case IV examines the effect of shadows cast by plants on soil. The solar 
and view angles are arbitrarily assumed to be 30* and 0* , respectively, and the 
reflectances are constant. The reflectance equation for a given pixel is 

R(A) = mR m ^(A) + gjRg^A) + ggR (A) (5.12) 

The scattergrams associated with the three aggregation levels for the Case 
IV simulation are shown in Figures 5.9-a,b,c. They reveal several interesting 
relations among percent cover and shadow, the level of aggregation, and the 
characteristic shape of the scattergram. 

All the data pairs fall within a space defined by a triangle. This is 
illustrated using the level 5 aggregation as indicated in Figure 5.10 (an expanded 
version of 5.9a). The vertices of the triangle (labeled Points B, C, and D) 
correspond to the assumed pure spectra of the full shadow (reflectance = 0.0), full 
canopy, and pure soil, respectively, as indicated in Table 5.1. 

Since pixels located within the triangle are linear mixtures of the three cover 
types, the exact percentage of any cover type can be determined on the basis of 
its location in the scattergram. For instance, the percent covers for an arbitrary 
pixel A shown on Figure 5.10 can be determined graphically as follows. First, 
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RED REFLECTANCE (%) 


Figure 5.10 Interpretation of scattergram, Case IV simulation, level 5 

aggregation: variable percent cover, constant soil and vegetation 
reflectances; shadowed soil background. 
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lines EF and GH are drawn through pixel A parallel to CD and BD, respectively. 
It is then noted that line EF is located about one-third of the distance between 
the line CD and point B. That indicates that pixel A contains 33 percent shadow. 
Line GH is situated about one-fourth the distance between the line BD and point 
C, indicating that pixel A has 25 percent vegetation cover. The remaining cover, 
42 percent, is bare soil, which can be checked on the basis of pixel A's location 
between line BC and point D. 

The above determination of the three cover types is simply a graphical 
illustration of an analytical solution applicable within the limits of the Case IV 
assumptions. It can be applied to any level of aggregation. The solution could 
also be achieved algebraically using equation (5.12) for both wavelengths (two 
equations) and equation (5.7). 

The Case IV scattergrams also reveal an important relation between shadow 
length scale and pixel size. For instance, at the level 5 aggregation, since the 
length scale of the shadow is about the same as the pixel scale, there are 
numerous instances when the shadow of a tree in one pixel falls onto an adjacent 
pixel. The three components of the pixel (vegetation, shadowed soil and 
illuminated soil) are independent of each other in a majority of cases. As a result, 
pixels can occupy almost any space within the limits of the triangular scattergram 
given a large enough sample size. 

As the level of aggregation increases, however, the length scale of the 
shadows becomes much smaller than the size of the pixel. As a result, shadows 
associated with a given tree fall increasingly within the same pixel and the amount 
of ground shadow becomes more and more correlated with the amount of 
vegetation cover. Mathematically, a covariance is generated among the three 
cover variables for the higher levels of aggregation which can be expressed, 
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(5.13) 

(5.14) 


h = 8 S ( m ) 

g, = g,(“) = 1 - m - g^m) 


and equation (5.12) becomes 


R(A) = R[A, m, gj(m), g g (m)] (5.15) 

A major consequence of the above relations is that it reduces the feasible 
region in the scattergram. Even at the level 5 aggregation (Figure 5.9a), that 
effect is manifested as a slight indentation in the upper right hand side of the 
triangular scattergram. At higher levels of aggregation, Equation (5.15) implies 
that there is only one position in the scattergram associated with a given canopy 
cover. As a result, one should expect the triangular scattergram observed at the 
level 5 aggregation to collapse to a single curved line when the shadow length 
scale becomes small relative to the pixel size. That is indeed shown to be true in 
a progressive manner by examining the sequential shapes of the scattergrams in 
Figure 5.9b (level 10 aggregation) and Figure 5.9c (level 30 aggregation). 

5.5 Case V - Shadowed S oil Background. Variable Soil Reflectance 

Case V is a more realistic version of the shadow model in which soil 
reflectance is assumed normally distributed as in Case II. The governing equation 
for an individual pixel is 


96 




The resulting scattergrams for the different levels of aggregation are shown in 
Figure 5.11. 

The scattergrams of the Case V simulation represent a combination of the 
effects illustrated in Case II (constant vegetation reflectance, variable soil 
reflectance) and Case IV (shadow effects). 

For instance, the scattergram of the level 5 aggregation, Figure 5.11a, 
exhibits a triangular shape overall, but with a pronounced indentation in the 
upper right portion due to the shadow effects. It can be regarded as a 
superposition of many triangular scattergrams, each for a homogeneous soil 
(constant background reflectance), similar to that of Case IV, level 5 aggregation 
(Figure 5.10). That is illustrated in Figure 5.12 (expanded version of 5.11a). 
Those triangles share two common vertices at 1) the point of pure shadow 
reflectance (point B), and 2) the point of pure vegetation reflectance (point C). 
The third vertex (labeled Dj, D2, D3, ... etc.) is unique for each triangle, 
representing the reflectivity of a particular soil which is homogeneous at that 
aggregation. The collection of all vertices, D, constitutes the true soil line. 

In the particular example shown, the true soil line has an intercept greater 
than zero, and is thus situated slightly inside the boundaries of the overall 
scattergram, as indicated in Figure 5.12. It is also possible, however, that the 
shadowed soil reflectance lies above the soil line. Only in such cases will the 
bottom of the scattergram accurately represent the true soil line. 

An important consequence of the level 5 aggregation is that pixels containing 
different mixtures of vegetation, shadow, and variable soil can occupy the same 


97 




98 


Figure 5.11-a,b,c Red-infrared scattergrams, Case V simulation: variable 

percent cover and soil reflectance, constant vegetation 
reflectance; shadowed soil background. 
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Figure 5.12 


Interpretation of scattergram, Case V simulation, level 5 
aggregation: variable percent cover and soil reflectances, constant 
vegetation reflectance; shadowed soil background. 



location in the scattergram. As a result, the percent cover of individual pixels can 
not be determined explicitly as shown in previous examples. 

The scattergrams of the levels 10 and 30 aggregation are shown in Figures 
5.11b and 5.11c, respectively. As in Case IV, because of the unique relation 
between shadow and vegetation cover at this scale, the scattergrams collapse 
progressively to the shape of a "tasseled cap" (Kauth and Thomas, 1976). At the 
level 30 aggregation, the scattergram consists of a series of juxtaposed curved 
lines, each line possessing constant average soil reflectivity (similar to Case IV, 
level 30 aggregation, Figure 5.9c), extending from individual points on the true 
soil line to the tip of the tasseled cap. That is illustrated in Figure 5.13 
(expanded version of 5.11c). 

Unlike the level 5 aggregation, percent cover can be estimated for Case V, 
level 30 aggregation, in a manner similar to Cases II and IV. Percent cover is 
proportional to the distance between the soil line and the tip of the tasseled cap. 
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Figure 5.13 Interpretation of scattergram, Case V simulation, level 30 
6 aggregation: variable percent cover and soil reflectance, constant 

vegetation reflectance; shadowed soil background. 
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Chapter 6 


REFLECTANCE AND COVER MOMENTS 


This chapter describes various moment equations applicable to the solution 
of the inverse problem described in Chapter 7. The moments include the mean, 
variance, cross-spectral covariance, and the spatial covariance of the reflectance 
equation. Conditional moments are formulated for portions of the scene where the 
pixels are assumed to possess one or more similar attributes, and which can be 
identified through ones knowledge of the scattergram. The geometric similarity 
and scaling criterion necessary for the application of the conditional moment 
equations are developed. A sampling scale ratio is developed as a quantitative 
scaling criterion to test when the fractional covers are functionally related. 



When one or more of the terms of the reflectance model are considered 
random, then the moments of both (4.4) and (4.5) can be expressed in terms of 
the moments of the individual variates. Such expansions can be achieved by 
applying fundamental properties of random functions without prescribing the 
probability density functions of the variates. For instance, if one assumes i) all 
the terms of equation (4.4) to be random, ii) that the reflectance terms are 
statistically independent of the percent cover terms, reflectances of different cover 
types are statistically independent, a general set of moment equations as presented 
below can be obtained. The mathematical details and the assumptions underlying 
the moment equations are provided in Appendix C. 
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( 6 . 1 ) 


The expected value, or mean pixel reflectance is 

E[R(A,x)] = £ E[f.(x))E[R.(A,x)] 
i 1 

where E[ ] designates the expected value. The variance of equation (4.4) can be 
written, 

VAR[R(A,x))j = E |E[fj(x)] 2 VAR[Rj(A,x)] + E[Rj(A,x)] 2 VAR[fj(x)] 

+ VAR[fj(x)] VARfRjC A,x)j | 

+ S E E[Rj(A,x)] EfR^A^jeov^Qc)^.^)] (6.2) 

1 t j 

where COV[fj(x),T(^)] represents the covariance between cover type i and j in a 
given pixel. The summations in (6.1) and (6.2) occur over the four cover types. 
The cross spectral covariance of the reflectance between the two bands, A, and A 
for a given pixel is written, 

COV[R(A 1 ,s),R(A 2 ,x)] = E {E[f.(x)] 2 COV[R i (A 1 ,x),R i (A 2 ,x)] 

+ VARfyx)] EfRjCApX)] E[Rj(A 2 ,x)]} 

+ E E E[R i (A 1 ,x)]E[R.(A 2 ,x)] COV[f.(x),f.(x)] 

* t j J 1 J 

(6.3) 

Similarly, the spatial covariance of the total reflectance between any two pixels in 
one band is, 
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COV x [R(A,x),R(A,x')] = S £ { [cOV x [f i (x),f:(x')] + E[fj(x)] E[f.(x')]' 

i j ll J i J 


COV x [R i (A,x),R j (A,x')] + E[Rj(A,x)] E[Rj(A,x')] 

- E[fj(x)] E[fj(x')] E[Rj(A,x)] E[Rj(A,x')] 

(6.4) 

where x = (xpX 2 ) an ^ x' = (x{,x^) represent two separate locations in the image 
and COV [R.(A,x),R.(A,x')] equals zero (by assumption) except for i = j. 

X 1 J 

The relationship among the moments of the cover variates can be established 
using equation (4.5). Those quantities are not assumed to be statistically 
independent but are related by the geometry and spatial distribution of the plants. 
The expected value is, 


l E[fj(jc)] = 1 (6.5) 

and the variance is, 

£ VAR[fj(x)] + E £ COV[f.(x),f.(x)] = 0 (6.6) 

i i # j J 

Equations (6.1) through (6.6) constitute at least six moment equations 
which, theoretically, can be augmented if more than one wavelength is used. For 
instance, equations (6.1), (6.2) and (6.4) each represent two equations when 
written for both the red and infrared bands. The actual number of moment 
equations depends on several factors, including the nature of the scene, the number 
of random variables in the model, and the linear independence of the moment 
equations at different wavelengths. 
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6.2 Conditional Reflectance and Cover Moments 


In some instances, conditional moments can be written for a portion of the 
pixels that possesses similar attributes. Examples include pixels possessing only 
bare soil, equal amounts of vegetation cover, or the same probability distribution 
describing plant spacing. The formulation of conditional moments for those cases 
generally reduces the complexity of the analysis provided that the appropriate set 
of pixels can be identified. 

One approach for identifying a set of pixels with common attributes is 
through the interpretation of multidimensional scattergrams. The previous chapter 
demonstrated, through the use of simulated images, that the structure of 
semi vegetated scenes manifests itself in the structure of red— infrared scattergrams. 
That knowledge of the structure of scattergrams provides a mechanism for 
identifying sets of pixels, not necessarily located within the same segment of the 
scene, for which conditional moments can be formulated. 

6 . 2.1 Soil Line Conditional Moments 

A relatively simple example of conditional moments is for the case of bare 
soil pixels. For many semivegetated scenes, bare soil pixels orient themselves 
along a preferred "soil line" at the base of triangular red-infrared scattergrams as 
described in Section 3 . 4 , or, 

R g/V *) = “ Kg/'W*' + 7 (6- 7 ) 

where a is the slope and 7 the intercept. In terms of equation ( 5 . 7 ), the only 
fractional cover type is bare soil, or 


g I “ 


1 


(6.8) 
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By conditioning the mean and variance of the reflectance equations (6.1 and 6.2) 
along the soil line, then, 


EpKAjdlg, = 1] 

= E[R (A*)] 

(6.9) 

var[R(A,3c) |g, = 1] 

= VAR[R (A,x)] 
6 I 

(6.10) 


resulting in two additional conditional moment equations for a given wavelength. 
The soil reflectance moments are related by, 


E[R gi (A m ,x)J = a E(R g] (A RED ,*)) + 7 

(6.11) 

VAR [Rgj(^ I R,-)] = Q E^gj^RED’-^ 

(6.12) 


6.2.2 Cover Moments for Statistically Homogeneous Spatial Distributions 

In many natural semivegetated regions, the spatial distribution of vegetation 
follows particular patterns which can be analytically prescribed. The statistical 
analysis of such spatial patterns in botany, with emphasis on coniferous 
vegetation, is well documented (Diggle, 1983). Knowledge of those distributions 
can be useful for relating the cover moments to the spatial distribution and 
geometry of the vegetation elements. 

One example is when a portion of the scene contains plants which can be 
assumed to follow a Poisson distribution in space. The mean of the fraction of 
illuminated vegetation for that segment of the scene can be written, 

E[m | Poisson] = 1 - exp[-/>(A t + A g )] (6.13) 
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where p is the Poisson spatial density, A fc is the average area of one plant canopy, 
Ag is the average area of the shadow cast by one canopy on the soil. The 
variance of the fraction of illuminated vegetation can be obtained in terms of m 
from the coverage problem derived empirically by Garwood (1947), or 


VAR[m | Poisson] 


m 


{1 


- 2AnT At/Ap 


->] 


+ 


2.3m 

k 


0.5 


(6.14) 


where A t is the average area of the plant, A is the area of the pixel, and 
15 ** 

k = (A t /A p ) ' . A graphical solution of (6.14) is provided in Figure 7.3 for the 
Case V simulation , level 10 aggregation, presented in Chapter 5. 


6.2.3 Geometric Similarity 

Many semivegetated landscapes possess only a few dominant species whose 
shapes can be represented by simple geometric figures, such as cylinders, cones, or 
spheres. The plants can be of different heights or sizes reflecting different stages 
of growth. In order to parameterize such shapes with a minimum number of 
variables, it is useful to assume that they are geometrically similar. Geometric 
similarity, employed in that sense, implies that the ratio of the plant height to 
some canopy width scale is a constant regardless of the size of the plant. In the 
case of conifer trees represented by cones, Li and Strahler (1985) assumed that the 
apex angle was constant. That assumption can be generalized to other geometric 
shapes as well. For instance, the similarity parameter of cylinders is the aspect 
ratio, b, defined as the ratio of the mean width, D, to mean height, H. For 
spheres, no similarity parameter is required, as the ratio of height to width is 
unity. 
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The geometric similarity assumption can be extended to include the solar 
angle, and thus, becomes particularly useful for establishing the analytical 
relationship between canopy area and shadowed area for different spatial 
arrangements. When no overlapping effects are considered, the shape of the 
plant's shadow is dictated precisely by the geometry of the plant itself and the 
solar zenith angle, 9. By defining rj as the ratio of shadowed area, A g , to plant 
area, A fc , or, 

V = Ag/A t (6.15) 


then for example, for square cylinders, 

tan0 

’ = T> 


( 6 . 16 ) 


Similarity parameters and corresponding rfs for different geometric shapes are 
provided in Table 6.1. A graphical illustration of geometric similarity in the case 
of cones is presented in Figure 6.1. 

The practical advantage of defining r? is that it allows one to absorb all the 
geometric factors which relate canopy area to shadowed area into only one 
variable. Consequently, the landscape can often be parameterized without the 
limitation of having to specify cones, cylinders or another geometric shape. 

As the vegetation density or the solar zenith angle increases, the shadow cast 
by one plant can extend far enough to be overlapped by the canopy of an adjacent 
plant. In such cases the amount of shadow is a function of the spatial distribution 
of the plants as well as of their geometry. Shadowing can occur when the plants 
are arranged in homogeneous deterministic spatial distributions, such as for row 
crops or orchards, or in stochastic distributions, as for natural vegetation. 
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Table 6.1 

Similarity of Canopy Geometry 


Canopy 

Shape 

Geometric 

Similarity 

Parameter 

Solar-Geometric 
Similarity Parameter 

i = V A t 

Circular 

Cylinders 

b = D/H 

4 tan0 
i 6~ 

Square 

Cylinders 

b = D/H 

tan0 

Cone 

<t> = tan _1 (D/H) 

(cot x ~ j + x )/ 7r 

Sphere 

none 

tan0 sin0 


D = mean canopy width 
H = mean canopy height 
9 = solar zenith angle 
X = sin'^tan 0/tan 0) 




For random spatial distributions, the relationship among cover types is more 
conveniently expressed in terms of their expected values. For example, in the 
particular case of Poisson distributed plants, equations (6.5), (6.13) and (6.15) can 
be combined to yield, 

E[gj] = jl - E[m]J^ + 1 (6.17) 


The derivation of equation (6.17) is provided in Appendix C. 

6-2-4 Conditional Moments for Pi x els with Constant Vegetation Cover 

The soil line moment equations are, in fact, a special case of the conditional 
moments for pixels of constant vegetation cover. For example, it was empirically 
shown in Section 5.4 that when i) the only variables in the scene are the 
fractional covers (vegetation, shadow, and illuminated soil) and soil background 
reflectance, and ii) a unique functional relationship exists among the different 
cover types of the form 

g I = 6i( m ’Ss) (6.18) 

then all pixels falling on a line parallel to the soil line possess equal amounts of 
vegetation cover. However, the distance of that line from the soil line was not 
linearly proportional to the amount of vegetation, but depended on the amount of 
vegetation and shadow, and the magnitude of the reflectances. 

Along each parallel line, since m, gj, and gg are constant, the conditional 
mean and variance of the reflectance equations given in (6.1) and (6.2) become, 
respectively, 


ill 


E[R(A,x) |m] = mR m (A) + g g R g (A) + gjE[R (A,x)] 

S I 

VAR[R(A,x) | m] = gj VAR[R (A,x)] 

6 I 


(6.19) 


( 6 . 20 ) 


The above formulation does not include the possible covariance between the soil 
reflectance and the amount of vegetation cover. Realistically, vegetation detritus 
changes soil reflectance by altering its organic and moisture content, and the 
resulting covariance must be considered in a more detailed analysis. 

6.2.5 Sampling Scale Ratio 

It is useful to develop a quantitative scaling criterion to test when the 
fractional covers are functionally related as in (6.18). One approach is to examine 
the relative scales of the shadow, as determined by the plant geometry and solar 
angle, and that of the pixel as determined by the field of view and altitude of the 
sensor. For instance, it can be reasoned based on the Case IV simulations, that as 
the scale of aggregation increases relative to the scale of the shadow, then, for 
homogeneous regions, a correlation develops among the different cover types. 

That correlation occurs since the variance in shadow cover is inversely 
proportional to pixel size. For very large pixels, the variance becomes so small 
that the fractional covers of each pixel approach the functional relationship in 
(6.18). The exact relationship depends on the geometry and spatial distribution of 
the plants, the solar angle, and the sensor characteristics. 

The above reasoning can be examined quantitatively by comparing order of 
magnitude estimates of the standard deviation (square root of variance) in shadow 
and the amount of canopy cover as a function of canopy geometry and pixel size. 
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For instance, in order for (6.18) to exist, one criterion that must be true is 

E[m] » VAR[g s ] 1/2 ( 6 .2i) 

From (6.6), (6.14), and (6.15) it can be shown empirically that 

VA R[gg] 1/2 ~ ^ (6 .22) 

P 

Thus, by combining (6.21) and (6.22), and noting that m is of the order of 
magnitude 10' 1 , a sampling scale ratio for Poisson distributions, S p , can be 
obtained of the form 



For Sp >> 10, equation (6.18) is valid. When vegetation reflectances are also 
constant, then (6.19) and (6.20) can serve as approximations for the more complex 
expressions (6.1) and (6.2). For homogeneous Poisson distributions, large S p 
implies that the partition of fractional covers in any given pixel approaches the 
mean relationship in (6.17), or 

gj s (1 - m) 7?+1 (6.24) 

For regular geometric (non— statistical) spacings, variability in g- occurs 
when the shadow associated with a given plant falls on a different pixel than that 
in which the plant is located. That situation is likely to occur when the scale of 
the plant is about the same or greater than the scale of the pixel. One criterion 
that avoids that situation is, 


A p >> ;?A t 


(6.25) 
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The sampling scale ratio for regular geometries is thus, 


S G “ lj$r t » 1 (6.26) 

The importance of (6.23), (6.24) and (6.26) for the inverse problem will be 
demonstrated in Chapters 7 and 8. 


63 Moments of Satellite Observed Radian™. 


It is useful to write the moments of the coupled landsurface-atmosphere 
radiation equations provided in Section 4.2 in order to examine the influence of 
the diffuse scattering terms on the moment equations. For instance, assuming 
horizontally homogeneous atmospheric conditions, the expected value of the 
satellite observed radiance given in equation (4.12) is 


E[L(A,x)] 



(6.27) 

where the L' values, defined in Section 4.2, are treated as constants for a given 
wavelength. The variance of (4.12) can be written 

VAR[L(A,s)]= 4 2 (A)VAR[Ef i (*)R j (A,i)] + l; 2 (A)var[s f!(i)a j (A,ic)] 

+ 2L;(A)L;(A)cov[[s f i ( J£ )R i (A,x)],[i; 


(6.28) 
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where, for simplicity, the summed terms have not been expanded as in (6.2) and 
(6.3). Equations (6.27) and (6.28) indicate that the atmospheric scattered diffuse 
radiance, L ( j(A), adds a constant term to the mean equation and has no influence 
on the variance. However, the diffuse radiance caused by adjacency effects can be 
significant, especially with regard to the variance, due to the introduction of 
cross-covariance terms. When the assumption that the surrounding area and 
target pixel reflectances are approximately equal, as given in (4.13), then the 
moments become 

E[L(A,x)j = [e E[fj(2t)]E[Rj(Ajc)]] L^(A) + L fl (A) 

VAR[l(A,i)] = L^(A)VAr(e fj(s)R j (A,2£)) 

In terms of the DN values recorded by the satellite, 

E[DN(A,x)] = E E[fj(x)]E[R|(A,x)] 4- l^A) (6.31) 

VAR[DN(A,x)] = VAR E fj(x)R!(A,x)| (6.32) 

where Rj(A,x) is the effective component reflectance defined in Section 4.3. Thus, 
it is shown that for regions of low interpixel contrast, the moment equations for 
the observed radiance, (6.31) and (6.32), are very similar in form to those of the 
target reflectance defined in Sections 6.1 and 6.2. That conclusion is valid only 
within the assumptions of the radiance model for optically thin atmospheres 
provided earlier. 


(6.29) 

(6.30) 
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Chapter 7 


GENERAL INVERSE PROBLEM: 

ESTIMATION OF SUBPIXEL PARAMETERS 

In most remote sensing applications, only the spatially integrated 
multispectral observations are available, with limited knowledge concerning the 
physical structure of the scene. That is especially true in natural areas where 
ground truth (e.g. training samples or spectral signatures) is not regularly 
obtained. 

The following sections present an approach for estimating the bulk physical 
parameters of the scene, with emphasis on subpixel vegetation cover, using the 
method of moments in the red and infrared bands. The method consists of 
equating the theoretical moments derived above for the two bands to the sample 
moments of the red and infrared images, and solving for the unknown parameters. 

Several versions of the method are possible depending on the nature of the 
scene. In this chapter, the method is applied to two simulation cases presented in 
Chapter 5. Inversion is first applied to Case II as an introduction to the method 
on a relatively simple scene, without shadows, in which the Sampling Scale ratio 
for Poisson distributions, Sp, is zero. The inverse procedure is then applied to 
two different aggregations of the same landscape simulated in Case V. Different 
versions of the inverse procedure are used on each aggregation of Case V, as a 
direct consequence of the change in value of the Sampling Scale Ratio with pixel 
size. 

The simulated images are analyzed primarily as a theoretical demonstration 
of the method for scenes in which the surface parameters are well known and can 
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be controlled. The analysis shows to what extent subpixel parameters, evidently 
lost through the aggregation process, can be retrieved. 


7.1 Inversion of Case II 

The parameters of the Case II image can be retrieved by employing the 
conditional moment equations to sets of pixels which possess similar attributes, as 
identified in the scattergram. More than one approach is possible, depending on 
which parameters are desired. This section demonstrates a version applicable for 
estimating the fractional vegetation cover on a pixel by pixel basis. The approach 
uses the knowledge that the sampling scale ratio, S p , is very large, since in the 
case of no shadows, equation (6.15) yields t ) = 0 and equation (6.23) yields 
S p = oo. It also takes advantage of the knowledge that pixels of constant 
vegetation cover lie parallel to the soil line. Once a set of pixels has been 
identified, the solution can be obtained by writing conditional moments along 
those lines using either band. 

The solution procedure is first, to locate the soil line in the scattergram, and 
then to calculate the sample mean and variance of those pixels. Equating those 
sample moments to the theoretical conditional moments of the soil line, equations 
(6.9) and (6.10), provides a direct estimate of the mean and variance of the soil 
reflectance at that scale. Next, a narrow band of pixels lying at an arbitrary 
distance from the soil line, but parallel to the soil line, is chosen and the sample 
moments of those pixels are calculated. The fractional vegetation cover of pixels 
in that band is estimated by equating the sample variances to the theoretical 
variances, and using (6.20) rewritten for Case II in which m = 1 - g , or 


m = 1 - 


VAR[R(A,x) | m]/VAR[R (A,2 c)]1 1/2 

J 


(7.1) 
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By choosing another parallel line at a different distance from the soil line, another 
estimate of m can be obtained for a new set of pixels. The procedure is repeated 
until the desired number of pixels has been analyzed. The computer code for the 
solution of the inverse problem for Case II is provided in Appendix D. 

The vegetation reflectance can be retrieved although it is not a prerequisite 
to estimate m. It is obtained by rewriting (6.19) specifically for Case II in which 
gs = 0, or 

R m (A) = (E[R(A,x) |m] - E[R (A,x)]}/m + E[R (A,x)] (7.2) 

©I 

The results of the Case II analysis using either the red or infrared band are 
shown in Figure 7.1 and Table 7.1. Figure 7.1 contains a plot of the estimated 
values of m versus the simulated values of each pixel in seven different parallel 
lines. The results indicate excellent retrieval for most values of m, with a 
standard deviation of error, s, equal to 0.026 for calculations using the red band 
and 0.028 for calculations using the IR band. The estimated reflectances also 
agree closely with the actual values for all values of m as can be seen in Table 7.1 
for one arbitrary line at which m equals 40 percent. 

It is noted that the above estimates of vegetation cover are made without 
introducing any assumptions on the geometry or spatial distribution of the trees. 

It is further noted that the variance of the soil reflectance computed above 
represents the variance of the aggregated process. The variance of the point 
process can be retrieved by applying the appropriate variance function to the 
results obtained at the aggregated level (Vanmarcke, 1983). 

For homogeneous regions, the spatial correlation function of soil reflectance 

can also be retrieved using the covariance equation (6.4) which, in Case II, reduces 
to, 
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ESTIMATED SUBPIXEL CANOPY COVER 



SIMULATED SUBPIXEL CANOPY COVER 


Figure 7.1 Estimated versus simulated canopy cover 
Case II simulation. 
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Table 7.1 

Comparison of Simulated and Estimated Scene Parameters 
Cases II and V, m = 0.40 

CASE II CASE V 




Simulated 

Estimated 

Simulated 

Estimated Values 

Parameter 

Units 

Value 

Value 

Value 

Method 1 

Method 

Reflectance Parameters 







R J X red) 

% 

15.0 

15.7 

15.0 

22.2 

14.8 

W 

% 

40.0 

40.6 

40.0 

51.3 

39.8 

R g s ( \ED ) 

% 



0.0 

o.of 

o.of 

\ |A ir> 

% 



0.0 

o.of 

o.of 

E t R g/ A RED» 

% 

15.0 

15.0 

15.0 

15.3 

15.3 

ElyV] 

% 

20.0 

20.0 

20.0 

20.3 

20.3 

VAR[R (A red )], level 10 

% 2 

4.8 

5.3 

4.8 

3.3 


VAR[R (A,^.,,)], level 30 

% 2 


.... 

3.0 


3.0 

VAR|R s (Aj R )), level 10 

% 2 

4.8 

5.3 

4.8 

3.3 


VAR[R s (A [R )|, level 30 

% 2 


.... 

3.0 


3.0 

Geometric Parameters 







Canopy Area, A^ 

2 

m 

1.0 


1.0 

0.85 


Similarity Parameter, 7) 


0.0 


2.0 

2.2 

2.0 

Sampling Scale Ratio, S 



00 


50 

450 

Soil Line Parameters 







Spatial Correlation 







Length Scale 

m 

20.0 

21.0 




Soil Line Equation 

% 


Y a ir> = 10 R 

g I (A RED ) 

+ 5.0 



T 


Assumed 


value 
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(7.3) 


COV x (R(A,i),R(A,js')] = E(g ] 2 + VAR(m] COV [R (A,x),R (A,s')l 

L J g I g I 

+ VAR[m][R m (A)-E[R gi (A,x)]] 2 

Two steps are required. First, the sample spatial covariance of soil 

reflectance, COV [R (A,x),R (A,x')], must be determined by computing (7.3) at 
Si g! 

different lags, x - x'. The sample spatial correlation function is then obtained 
from 



to which an appropriate function can be fitted if an analytical relationship is 
desired. 

The correlation function for Case II was estimated by applying the above 
procedure to the eight different segments (representing eight values of m) of the 
red and infrared scenes. The results are shown in Figure 7.2 and Table 7.2. 

Figure 7.2 shows a comparison of i) the sample correlation of the simulated bare 
soil reflectance in the red band at the 10-meter level of resolution and ii) the 
estimated soil reflectance correlations using (7.3) and (7.4) for a range of m's (14, 
40, and 57%) at lags of 10, 20, 30 and 40 meters. The good agreement between 
the simulated and estimated correlations, as indicated in Figure 7.2, was typical 
for all values of m in the red band, but poor for the IR band for m greater than 
about 50%. 

The estimated soil reflectance length scales for each value of m, assuming an 
exponential correlation, are shown in Table 7.2 for both the red and infrared 
bands. Overall, the estimated values compare favorably with the actual simulated 
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LAG (METERS) 

Figure 7.2 Estimated soil reflectance spatial correlation compared to sample 
soil reflectance correlation, Ca$e II, red band, level 10 aggregation. 




Table 7.2 


Estimated Soil Reflectance Length Scale 
Level 10 Aggregation 

Simulated Length Scale = 20 m 


Fractional 

Canopy 

Cover 

IR 

Case II 

RED 

IR 

Case V 

RED 


(m) 

(m) 

(m) 

(m) 

0.05 

20.6 

21.4 

21.2 

20.7 

0.14 

16.8 

21.3 

18.4 

20.5 

0.26 

22.5 

21.0 

12.3 

13.7 

0.39 

18.4 

21.0 

— 

14.9 

0.51 

20.0 

21.3 

— 

— 

0.52 

12.2 

21.2 

— 

— 

0.63 

11.3 

21.07 

— 

— 

0.78 

18.3 

20.3 

— 
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value (20 m), although the agreement is better for the red band. In the infrared 
band, the length scales are not retrieved as well at higher values of m. That is 
due to the greater difference in the reflectances of the vegetation and soil in the 
IR band, as compared to the red band. 

7.2 Inversion of Case V 

For the Case V simulation, equation (6.16) yields rj = 2.0. Two different 
approaches are presented for estimating the Case V parameters, depending on the 
magnitude of the sampling scale ratio, S , defined in (6.23). A first approach 
(Method 1), applicable for all ranges of Sp, consists of solving the full set of 
moment equations simultaneously for different statistically homogeneous regions. 
The second procedure (Method 2), applicable only for large Sp in which (6.18) is 
true, utilizes a simpler approximate set of moment equations and one's knowledge 
of the structure of the scattergram. 

7.2.1 Estimation of Parameters. Method 1: Sp 10 

The introduction of shadows in Case V adds to the complexity of inverse 
problem in several ways. First, it increases the number of cover types to three 
(e.g., illuminated canopy, illuminated soil background, and shadowed soil 
background) as well as the number of associated reflectance terms in both bands. 
Additionally, the covariance among the three cover types must be considered in 
the analysis. 

At the level 10 aggregation, the sampling scale ratio, Sp, is only ^ 10 - 
Since (6.24) is not valid for such small Sp, pixels with different amounts of 
vegetation cover can occupy the same location in the scattergram, previously 
shown in Figure 5.9b. The assumption of constant canopy cover for pixels 
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oriented parallel to the soil line does not apply (except for the soil line itself in 
which m = 0). Thus, unlike Case II, subsets of pixels required for the analysis 
can not be identified on the basis of the scattergram. 

When segments of the scene can be assumed statistically homogeneous, 
however, the following conditional moments can be written for the set of pixels 
located in each region. The mean reflectance equation is, 

E[R<Ajt)J = E[m]R m (A) + E[gg]R g (A) + E[gj]E[R (A,fl) (7.5) 

S ®I 


The variance equation now includes the variances and covariances of the individual 
cover types, or 

VAR[R( A, 2 C) ] = E[gj] 2 VAR[R (A,x)] + R m (A) 2 VAR[m] + E[R (A,x)] 2 VAR[g ] 

©I ©| I 

+ VAR[ gl ]vAR[R gj (A,x)] + R gs (A) 2 VAR[g s ]-2R m (A)E[R g (A,x)] ♦ 

jvAR[m] + COV[m,g s ]} + 2R m (A)R g (A) COV[m,g s ] 

s 

- 2E[R gi (A, 2 c)]Rg s (A){vAR[g s ] + COV[m,gg] j (7.6) 


The cross covariance between any two spectral bands is written, 
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C° v (R(A red , s ),R(Air,x)J = VAR|m]R m (A RED )R m (A IR ) 

+ E[gj]COV[Rg^(A RED ,x),Rg^(A IR ,x)] 

+ VAR[g I ]E[Rg^(A RED ,x)]E[Rg^(A IR ,x)] + VAR[gg]R^(A RED )R^(A IR ) 

“ { R m^RED^ E f R gi ^lR’ x )] + R m^IR^ E t R g^ A RED’-^} 

• |vAR[m] + COV[m,g s ]} 

+ { R rn^RED^ R g s ^IR^ + R m^IR^ R gg^RED^} COV [ m ig s ^ 

“ { E ^ R g/^RED’ X ^ R g s ^IR^ + E f R gl ^IR’-^ R g s ^RED^} 

• {vAR[gj] + COV[m,g s ]} (7.7) 


The derivation of (7.5), (7.6), and (7.7) are given in Appendix C. The moments 
of the percent cover variates can be written, 

E[m] + E[gj] + E[g s ] = 1 (7.8) 

VAR[m] + VAR[g s ] - VARfej] + 2COV[m, gg] = 0 (7.9) 

for a total of seven conditional moment equations when (7.5) and (7.6) are written 
for two bands. The addition of the five soil line equations, i.e., (6.9) and (6.10) 
written for both bands, and (7.7) written for m = 0, brings the total to twelve. 
The unknowns include i) the means and variances of the three cover types plus 
COV[m,g s ], ii) the vegetation and shadowed reflectances written for both bands, 
iii) the mean and variance of the illuminated soil for both bands, and the illu- 
minated soil reflectance spectral cross covariance. The number of unknowns is 16. 
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The solution procedure is, first, as in Case II, to obtain the four illuminated 
soil reflectance moments using (6.9) and (6.10). The cross covariance of the soil 
reflectance is also directly obtained from those moments or, equivalently, from 
conditioning (7.7) along the soil line. 

Since there are still four more unknowns than equations, additional 
assumptions of i) Poisson spatial distribution and ii) geometric similarity of the 
plants are required. Those assumptions introduce two unknowns, 77 and A , but 
provide four additional equations, namely (6.14), (6.17), and 

VAR[ gj | Poisson] = h 2 {E[ gl ],^,A t ,A p } (7.10) 

and 

COv[m,gg] = h 3 {E[m],i,,A t ,A p } (7.11) 

Graphs of hj (equation 6.14), h,, and h 3 , are shown in Figure 7.3 for the Case V 
simulation, level 10 aggregation. 

The covariance between m and g g was obtained in the following manner. It 
was observed during the simulations that the correlation between m and g could 
be approximated by a stepwise linear function of m, as shown in Figure 7.4. The 
expression for COV(m,g g ) in (7.11) and Figure 7.3 was obtained by combining (7.9) 
and the definition of correlation, or 

COV(m,g ) 

^m,g ~ ‘ (7-12) 

[COV(g s )VAR( gl )]°- 5 

where 


p = 

1.00 - 0.50E[m] 

for 0.00 < E[m] < 0.10 

p = 

1.28 - 3.27E[m] 

for 0.10 < E[m] < 0.65 

p = 

-0.57 - 0.43E[m] 

for 0.65 < E[m] < 1.00 


(7.13) 
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FRACTIONAL CANOPY COVER, m 


Figure 7.4 Functional relationship 



versus m for Case V simulation. 
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in order to yield the following semiempirical expression 


COV 2 (m,g ) VAR 2 (g ) - VAR 2 (m) 

0 = 5 + cov(m,g ) 

2 p VAR (m) 2 

(7.14) 


which can be solved by iteration for COV(m,g g ). 

Since there are still two more unknowns than equations, the solution to the 
inverse problem can be obtained in one of two manners. One approach is to 
choose two different homogeneous regions, and to write conditional moment 
equations for those two regions. Each segment adds seven unknowns (the cover 
moments), but provides nine equations ((6.14), (6.17), and (7.5) through (7.9), with 
(7.5) and (7.6) written for two bands), for a net gain of two. Thus, theoretically, 
the parameters can be obtained by simultaneously solving 25 conditional moment 
equations simultaneously. 

A second, more practical approach in some cases is to assume that R (A) in 

g S 

equations (7.5), (7.6) and (7.7) is negligible compared to the other terms. That 
assumption is only valid for i) relatively clear atmospheric conditions when the 
diffuse radiation is small, ii) low zenith angles (i.e., small rj), and iii) bands in 
which the vegetation and illuminated soil reflectances are both much greater than 
the shadowed soil reflectance. The third condition is generally true for most red 
bands, and for near-infrared bands when the extinction due to leaf area is large. 
The solution can thus be obtained using one set of conditional moment equations 
applicable to one homogeneous region. 

Due to the complexity of the equations, the solution to the Case V, 

Method 1 inverse problem can not be solved explicitly. Instead, it is obtained by 
minimizing the sum of the squared errors between the theoretical and sample 
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moments. The procedure involves iteration over E[m], A., and 77 until the 
minimum error is obtained. 

Since the complete set of equations is non-linear, in order to avoid local 
minimum solutions, the global solution is obtained by solving the problem over a 
range of reasonable initial values and iteration steps. The computer code for the 
Case V, Method 1 inverse problem is provided in Appendix D. 

The results of applying the second approach to the level 10 aggregation of 
Case V are shown in Figure 7.5 and Table 7.1. The procedure provides estimates 
of only the mean values for each statistically homogeneous area. The results 
indicate that the values of the estimated mean vegetation cover compare 
reasonably well to the simulated mean values (the standard deviation of error, 

s = 0.056). The reflectances are also fairly well recovered, although A and A { are 
not. 

The length scales of the spatial correlation function of soil reflectance were 
also computed for Case V - Method 2 using (7.3) and (7.4) for various m's. The 
results, shown in Table 7.2 indicate that the retrieval of such length scales was 
only possible for low values of m. The poor retrieval is a result of the large 

amount of shadow which masks the soil background in this simulation at large 
m's. 

7 - 2 - 2 Estimation of Paramete rs. Method 2: S p jn 

At the level 30 aggregation, S p equals ^ » 10 . A t such large S p , an 
approximate functional relation occurs among g { g g and m as shown in (6.18). 
Since the simulated canopy reflectance is also constant, then pixels with different 
amounts of canopy cover will orient themselves parallel to the soil line as 
described in Sections 5.4 and 6.2.4, and shown in Figure 7.6. Thus, the inverse 
procedure for large S p makes use of that knowledge by formulating approximate 
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ESTIMATED MEAN SUBPIXEL CANOPY COVER 



SIMULATED MEAN SUBPIXEL CANOPY COVER 


Figure 7.5 Estimated mean versus simulated mean subpixel canopy cover, 
Case V simulation, level 10 aggregation. 
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INFRARED REFLECTANCE (%) 



RED REFLECTANCE (%) 

Figure 7.6 Indication of direction of increasing canopy cover, 
Case V simulation, level 30 aggregation (method 2) 



conditional moments for individual subsets of pixels which lie along those parallel 
lines. 


As in Case II, once the appropriate set of pixels needed for the conditional 

sample moments are determined from the scattergram, the parameter estimation 

can proceed using only one band. For each arbitrarily chosen parallel line, there 

are seven unknowns including i) the three fractional covers m, gj, and g g , 

ii) two constant reflectance terms, R (A), R (A), and iii) the mean and 

g s 

variance of the soil reflectance, E[R (A,x)], and VAR[R (A,x)], respectively. 

§1 gj 

The available conditional moment equations are i) the two soil line equations, 
(6.9) and (6.10), ii) the conditional mean reflectance, (6.19), and iii) the 
variance as provided in (6.20). It is noted that (6.20) is obtained directly from 
(7.6) by assuming that the variances of the fractional covers are negligible. 

In addition, the percent cover equation, (5.7), gives, 


m + g s + % = 1 (7.15) 

for a total of only five equations, two less than the number of unknowns. The 
above formulation only allows several terms to be retrieved. For instance, the soil 
line moments can be obtained as in the previous examples. Then, by selecting an 
arbitrary locus of pixels parallel to the soil line, gj can be solved directly using 
(6.20). However, close examination of (6.19) reveals that m cannot be determined 
using only the above equations. 

The remaining terms can be obtained by first assuming Poisson spatial 
distribution and geometric similarity of the plants. Those two assumptions permit 
one to use (6.24), for a total of six equations, although an additional unknown, 77, 
is introduced. Since there are still two more unknowns than equations, the deficit 
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is made up in one of two ways. One approach is to arbitrarily choose two 

additional parallel lines (for a total of three) for which the above conditional 

moments apply. Each additional line introduces three new unknowns (m, g , and 

s 

gj), but augments the number of equations by four (equations (6.19), (6.20), 

(6.24), and (7.15)). Thus, m for each line is obtained by first solving for gj for all 
three parallel lines using (6.20), and then solving the remaining moment equations 
simultaneously. While this approach is theoretically correct, its validity is limited 
to scenes which possess a large number of pixels in at least three different 
homogeneous regions. 

As in Method 1, a second approach is to assume that the shadowed 

reflectance term in (6.19), g_R (A), is negligible compared to the other two 

5 6 S 

illuminated terms. That assumption eliminates the need to calculate R (A) and 

«s 

the analysis can be conducted using only two parallel lines. 

The results of the second approach are summarized in Figures 7.7 and 7.8 
and Table 7.1. Figure 7.7 contains plots of the estimated canopy cover versus the 
simulated values for each pixel contained in each of five conditional lines, using 
both the red and infrared bands. The standard deviation of error s, is 0.028 in 
the red band and 0.069 in the infrared band. Although the agreement is very 
good, Figure 7.7 indicates that the estimated values are generally lower than the 
simulated values, especially at higher values of m. That difference is due to the 
error associated with neglecting the cover variance terms in (7.6). The error is 
greater for estimates made using the infrared band since the magnitudes of the 
infrared reflectances are greater than those in the red band. 

Figure 7.8 indicates equally good agreement between the simulated values of 
gj(m) and g g (m) and the theoretical curves using (6.24), (7.15) and the mean 
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estimated value of 77 . The simulated and estimated reflectance terms also compare 
favorably, as shown in Table 7.1. 
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ESTIMATED SUBPIXEL CANOPY COVER 



Figure 7.7 Estimated versus simulated subpixel canopy cover, 
Case V simulation, level 30 aggregation (method 2). 
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SUBPIXEL ILLUMINATED SOIL 
OR GROUND SHADOW 



Figure 7.8 Comparison of simulated values of illuminated and 
shadowed soil fractional cover versus canopy cover, 
to the theoretical curves using equations (6.24) and' 
(7.15), Case V simulation, level 30 aggregation. 
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Chapter 8 
CASE STUDIES 


This chapter examines the application of the canopy reflectance model to 
actual multispectral data obtained over two field sites. The first site is a pecan 
orchard in southern Arizona for which actual aerial radiometric data were 
obtained. The second site is a pinyon-juniper watershed in northern Arizona for 
which both aerial and satellite multispectral data (Landsat TM data) were 
obtained. Atmospheric effects on the subpixel estimates are examined using 
hypothetical values of backscattered solar diffuse radiation. 

8.1 Scattereram of a Pecan Orchard 

The pecan orchard represents a special case of the canopy reflectance model 
in which the trees are spatially distributed in a fixed geometric fashion and the 
only random property is the soil background reflectance. In this example, a visual 
comparison is made between the plots of the radiometric data in the red-infrared 
reflectance space, and a hypothetical scattergram constructed from ground truth 
measurements at the time of overflight. The moment analysis is not applied to 
the estimation of subpixel cover due to the limited number of pixels. Rather, a 
qualitative comparison of the theoretical and estimated canopy cover is made. 

8.1.1 Site Description 

The study site is located within a flat one mile square area near Maricopa, 
Arizona, about 40 km south of Phoenix. Aerial radiometric measurements were 
collected at an altitude of about 150 meters at 9:30 a.m. on June 12, 1988, as part 
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of the MAC (Maricopa Agricultural Center) III Experiment organized by the 
Water Conservation laboratory, Agricultural Research Service, Phoenix, Arizona. 
During the experiment, there were no clouds, and the air could be qualitatively 
described as clear and dry. 

The orchard itself consisted entirely of pecan trees planted on a square grid 
in an east-west orientation, with center intervals of approximately 85 meters. 

The diameter of individual trees ranged from about 5 to 10 meters (A ~ 20 - 80 

b 

sq.m.), with a height to depth ratio of about unity. Tree height was generally 
constant in any given section of the orchard, and thus, tree canopies were not 
significantly shadowed by adjacent trees. The size of the trees and, thus, the 
amount of canopy cover, could vary from pixel to pixel. Trees were interspersed 
with a combination of bare soil and senes ced grasses. 

8.1.2 Reflectance Data 

Radiometric observations were made using an Exotech radiometer with 
Thematic Mapper red (0.62-0.69 /nn) and infrared (0.78-0.90/m) filters at a 
ground resolution of about 40 meters (A p ~ 1250 sq.m.). The solar angle was 
estimated to be 43.5* at the time of overflight. 

Radiometric observations over the pecan orchard were converted to 
reflectance factors (ratio of target reflectance to the reflectance of a Lambertian 
surface; See Jackson et al, 1987) by Moran (1988). The reflectance factors are 
proportional to and approximately equal to the target reflectance, and thus for 
simplicity, they will be termed reflectance throughout this section. 

Ground truth values of the pecan canopy's bulk reflectance could not be 
easily obtained due to the large size of the trees. However, the aerial observations 
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over nearly continuous canopies indicate that the bulk canopy reflectance is about 
2 - 4% in the red band and 45 - 55% in the IR band. The mid-points of those 
values (canopy red reflectance = 3.0%; canopy IR reflectance = 50%) were 
arbitrarily chosen as estimates of the canopy reflectances. Shadowed reflectances 
were assumed to equal 10% of the canopy reflectances. 

Since no treeless pixels existed in the orchard itself, the soil reflectance and 
the soil line was obtained by sensing bare soil fields (Maricopa field No. 18, 27, 
lnd 32) immediately adjacent to the orchard which contained a mixture of bare 
soil and senesced grasses. The soil line obtained from a red-infrared plot of the 
data is shown in Figure 8.1. The line exhibits a nearly linear relationship as 
described in equation (3.5). The mean, standard deviation, and covariance length 
scale (computed as the average e-folding distance of the empirical correlation 
function) of those soil pixels, together with the parameters of the soil line are 
provided in Table 8.1. 

8-1-3 Fractional Cover Estimates 

As a substitute for ground truth, independent estimates of fractional cover 
were made by analyzing the histograms of the digitized multispectral video images 
for each of four radiometric observations. Because the length scale of the canopy 
is several meters, and the length scale of the video pixel is only about one-half 
meter, a majority of the pixels will be approximately pure canopy, pure shadow, 
or pure illuminated soil pixels. As a result, if the reflectance of each cover type is 
unique, a histogram of the digitized video image should possess local modes 
corresponding to each of the different cover types. The percentage of pixels 
associated with each mode approximates the amount of a particular fractional 
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Figure 8.1 Red-infrared scattergram of pecan orchard and bare soil. 
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Table 8.1 


Mean Reflectance (%) 

Standard Deviation 
of Reflectance {%) 

Soil Line Equation (%): 
Covariance Length Scale 


Soil Line Parameters 
of Pecan Orchard 


Red TM Band 
f0.62- 0.69 mu) 

27.3 

3.7 


Infrared TM Band 
(0.78- 0.90 mn) 

32.8 

4.1 


V A m 


) = 109 y w + 306 


200 meters 
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cover, depending on the wavelength. 

For example, Figure 8.2 contains the red and infrared histograms of one 
video image (time = 9.3969 hours) for which a radiometric observation was 
simultaneously made. In the red image, trees and shadows appear very dark, 
while the soil is bright. The corresponding red histogram exhibits a strong 
bimodal shape, with a local minimum occurring at an intensity level of 89. Thus, 
intensity levels less than 89 (31 %) are assumed to represent pixels containing 
primarily vegetation and shadow. Pixels greater than 89 are assumed to represent 
illuminated soil (69 %). 

For the infrared image, only the shadows appear dark while both trees and 
soil are bright. Since the amount of shadow is small, a strong bimodal effect is 
not observed in the histogram, although a slight trough is observed at an intensity 
level of about 128. As a result, pixels with intensity < 128 are assumed to 
represent shadow (12%), while pixels with intensity > 128 represent vegetation 
and illuminated soil (88%). Combining the results of the histogram analysis for 
both bands yields estimates of fractional cover for the vegetation (19%), shadow 
(12%), and illuminated soil (69%), for the single radiometric observation. 
Quantitative estimates obtained in the above manner were confirmed by visual 
examination of the video image. 

8.1.4 Comparison of Actual and Hypothetical Scattergrams 

The hypothetical scattergram is obtained by first conceptualizing the orchard 
as a stochastic geometric surface, in which the only variables are the fractional 
canopy cover and the soil reflectance. The hypothetical reflectance of a pixel was 
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assumed to be 


R(A) - mR m (A) + S s R gg ( A ) + SjRg^A,*) (8.1) 

where R (A), R (A) and R (A,x) represent estimated ground truth reflectances 
6 s 6 i 

of the illuminated canopy, shadowed soil, and illuminated soil, respectively. 

The next step is the calculation of the sampling scale ratio for regular 
geometries, Sq, based on the similarity parameter. Visual observations at the 
time of the experiment indicated that the trees could best be represented by 
circular cylinders, and that they were approximately geometrically similar. From 
Table 6.1, 


T) = 4H tan 0 (8.2) 

ttD 

Inserting the parameters of the experiment (0 = 43.6 degrees, H/D ~ 1) into the 
above equation yields rj = 1.21. The sampling scale ratio for geometric 
distributions is 


S G = V^t s 13 " 50 ( 8 - 3 ) 

Thus, based on the criterion given in (6.26) for geometric distributions, since 
Sq >> 1, the assumption of large Sampling Scale Ratio is made. 

The assumption of i) geometric similarity and ii) large sampling scale allows 
one to formulate a unique analytical relationship among the fractional covers as 
described in Section 6.2.4. For the particular orchard described above, two 
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different shadow regimes can be identified. Regime I occurs when the trees are 
small and the entire shadow cast by a tree is observed. In this case the fractional 
shadowed area g is linearly related to the fractional canopy cover, m or 

b 


for 



-1 < 0 


(8.4) 


where fi equals the tangent of the zenith angle and f is a similarity parameter 
equal to the ratio of canopy diameter to tree height. Regime II occurs for larger 
trees when the shadow cast by one tree extends far enough as to be overlapped, in 
part, by the canopy of an adjacent tree. A second term is added to the above 
equation to account for that decrease in shadowed area, or 



for 0 < 2 

m 

1 + 

- 1 < 2 


7T 

f 



In both regimes, the illuminated soil is constrained 

Si = 1 * m - h 


by 


(8.5) 


( 8 . 6 ) 
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A sketch of both regimes is provided in Figure 8.3 The graphical forms of 
g s (m) and g^m) are given in Figure 8.4. Also plotted on Figure 8.4 are the 
actual fractional cover estimates of several pixels obtained from aerial video. The 
plots indicate that the theoretical curves of the fractional shadow and illuminated 
soil agree reasonably well with the actual data. 

Using (8.1) through (8.6), the scattergram of a hypothetical orchard scene 
was constructed by superposing canopy cover ranging from 10 to 70 percent onto 
each of the soil background pixels. The resulting scattergram based on that model 
is shown in Figure 8.5. It possesses many similarities to the simulated cases 
presented earlier, including a triangular shape with curved sides and a flat base. 

Figure 8.5 also includes the plot of several radiometric data from Figure 8.1 
for which the subpixel fractional covers were estimated. The orchard itself does 
not possess a wide range of vegetation cover needed to establish a complete 
triangular scattergram as in the simulations. However, a comparison of the actual 
data with the hypothetical scattergram indicates that their location is consistent 
with the predicted values. A summary of the actual and hypothetical fractional 
covers for four pixels is provided in Table 8.2. The good agreement achieved 
above serves as a preliminary confirmation of the validity of the canopy 
reflectance model for explaining how subpixel variations in cover type affect the 
relative location of pixels in a red-infrared scattergram. 

8-2 Pinvon-Juniper Watershed: Aerial Radiometric Data 
This example tests the canopy reflectance model and inverse procedure on an 
actual semivegetated watershed for which aerial radiometric data were obtained. 


148 



Section 





Regime I Regime II 

(No overlap) (Canopy overlaps shadow) 


Figure 8.3 Hypothetical shadow regimes for pecan orchard. 
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SHADOWED OR ILLUMINATED SOIL (%) 



Figure 8.4 


Theoretical relationship of percent illuminated soil, g., and 
percent shadowed soil, g , as a function of percent canopy cover, 
m, compared to actual data for a pecan orchard. 




Figure 8.5 Hypothetical scattergram of pecan orchard compared with actual 

aerial radiometric data. 


151 



Table 8.2 


Comparison of Actual and Hypothetical 
Fractional Covers for 
Pecan Orchard 


Pixel 

Number 

Total 

Canopy 

Cover 


Total 
Shadowed 
Soil Cover 



Total 

Illuminated 
Soil Cover 

(%) 








Actual Model 

Actual 

Model 

Actual 

Model 

1 

55 

54 

35 

29 

10 

17 

2 

50 

51 

30 

29 

20 

20 

3 

15 

21 

11 

25 

74 

54 

4 

19 

20 

12 

24 

69 

56 
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It also shows how the procedure can be adapted to situations in which the total 
number of pixels is small. In such cases, the red-infrared scattergram may not 
form a fully developed tasseled cap shape, as previously shown for the idealized 
simulated cases. 

8.2.1 Site Description 

The study site is a small, natural semivegetated watershed, about 0.8 km 2 , 
located in the Beaver Creek Basin in the Coconino National Forest in north 
central Arizona, as shown in Figure 8.6. The area is relatively flat sloping 3.0 
percent to the southwest at an average elevation of 1900 meters. The 
predominant tree species is alligator juniper ( Juniperus deppeana ), a short, 
egg-shaped evergreen with tiny scale-like leaves, ranging in height from 3 to 5 
meters. Small amounts of Utah Juniper ( Juniperus osteosperma), a tree similar in 
shape to the Alligator Juniper, and Ponderosa Pine ( Pinus ponderosa), a taller, 
narrower evergreen with a rounded crown, also exist. The area between the trees 
is interspersed with a mixture of bare soil and a variety of sparse, relatively dry, 
semiarid grasses and shrubs. Field observations indicate that the fractional 
pinyon-juniper canopy cover ranges from about 0 to 70 percent, with a mean of 
about 25 percent. Soils are rocky and developed from volcanic materials, 
primarily basalts (Clary et al, 1974). 

8.2.2 Acquisition of Radiometric Data and Ground Truth 

Multispectral data were collected in a similar manner as for the pecan 

orchard experiment. The overflight occurred between 10:15 and 10:30 a.m. on 
June 23, 1988. at an altitude of about 150 m. using nadir-viewing instruments. 
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Figure 8.6 Location Map: Beaver Creek Watershed. 
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Approximately 200 radiometric observations were made using an Exotech 
radiometer with Thematic Mapper red (0.62-0.69 pm) and near infrared (0.78-0.90 
pm) filters at a ground resolution of about 40 meters. Simultaneous multispectral 
video was also obtained. 

Late June was chosen as the acquisition period as it was a relatively dry 
period when soil moisture was low and the grasses were in a somewhat senesced 
state, offering good contrast to the dark green evergreens. The solar zenith angle 
is relatively low compared to other seasons which minimizes the effect of shadows. 
Measurements were taken at 10:00 AM in order to conduct the analysis at the 
same time as a typical Landsat overpass, and to avoid further buildup of haze 
which was occurring during the morning of the overflight. 

The radiometer data, recorded in terms of voltages, were not converted to 
reflectances, as it would have required additional ground-based instrumentation, 
not typically available in most remote sensing applications. Further, since the 
primary interest was in estimating fractional cover amounts, conversion to 
reflectances was not necessary. Sensor voltage is approximately linearly related to 
incoming radiance (Jackson et al, 1987). Thus, in the analysis which follows, 
voltage was used as a surrogate measure of reflectance. The analysis is valid as 
long as i) the time interval over which the data were collected was small (several 
minutes) in order to minimize the effect of changing solar zenith angle, ii) the 
solar irradiance on all target pixels was constant, and iii) diffuse atmospheric 
effects on the reflected radiance were minimal and constant over the region. 

Estimates of fractional cover were made using the video data described in 
Section 8.1.3. The resulting histograms also exhibited a bimodal shape as for the 
pecan orchard, as shown in Figure 8.7 for one observation. 
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Figure 8.7 Histograms of red and infrared video images of the same 
semivegetated pinyon-juniper landscape. 



8-2.3 Estimation of Subpixel Canopy Cover 

As in the previous cases, the first step was to plot the entire set of 
observations in the red-infrared space as shown in Figure 8.8. The resulting 
scattergram, now in terms of voltages, possesses an overall triangular shape, 
although not as well defined as for the simulated cases. A relatively flat base does 
exist, however, the top of the scattergram is somewhat rounded and does not 
possess a fully developed "tasseled cap" shape. The reason for this lack of 
definition is that, unlike the simulated scenes, the watershed does not possess a 
full range of combinations of fractional vegetation cover and soil reflectance. 

The second step in the inversion procedure was to conceptualize the 
pinyon-juniper landscape as Poisson distributed spheres resting on a flat surface. 
The spherical trees are assumed to exhibit a constant bulk reflectance in each 
band, R m (A), which accommodates both illuminated and shadowed portions of the 
canopy. Soil reflectance is variable and can be shadowed or illuminated. The 
total reflectance of any pixel is thus given by 

R(A) = m R m (A) + g s R g (A) + gjR (A,x) (8.7) 

S I 


The geometric similarity assumption allows tj to be estimated directly from Table 
6.1, Item (iv). The value of rj was calculated to be 0.294 based on an estimated 
solar zenith angle of 30.25 degrees, obtained using Iqbal (Chapter 1, 1983). The 
mean tree diameter is on the order of several meters, thus, 

S P ~ 094A. >> 10 (8-8) 

L 
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INFRARED INTENSITY (VOL' 



RED INTENSITY (VOLTS) 


Figure 8.8 Red-infrared scattergram of pinyon-juniper watershed, 
aerial data. 
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and the assumption of a large sampling scale ratio is made. 

Next, the mean and variance of the soil background radiance (expressed in 
terms of voltage) were calculated by fitting a straight line through approximately 
20 pixels located at the bottom of the scattergram (Figure 8.9), and by estimating 
the moments of those pixels. The pixels do not fall completely on a straight line 
since they contain small amounts of green and senesced grasses, in addition to the 
bare soil. The estimated values of the soil moments are provided in Table 8.3. 

According to the method for Sp >> 10, the inversion procedure next 
requires the identification of sets of pixels lying in a band parallel to the soil line. 
The sample moments of those pixels are then used to estimate gj using (6.20). 
However, in the present example, the number of pixels in any given line is too 
small to compute sample moments, and a slightly different approach must be 
taken, requiring two steps. First, instead of choosing subsets of pixels lying 
parallel to the soil line, the entire ensemble of pixels (except those associated with 
the soil line) were analyzed simultaneously in order to obtain overall estimates of 
vegetation reflectance and fractional cover statistics for the watershed. The 
analysis required (6.17), (7.8) and approximate relations for the mean (7.5) and 
variance (7.6) based on the following order of magnitude analysis. 


Since tj = 0.294, then from (6.17) and (7.8) it is reasoned that E[g ] will be 

smaller than E[gJ and E[m] for m < 0.5. Since R (A) is likely to be less than 

% 

R (A) and E[R (A)], the product, E[g ]R (A), will also be small compared to 


g 


V 


the other two terms in (7.5). The mean reflectance can thus be approximated by 



(8.9) 
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INFRARED INTENSITY (VOLTS) 



Figure 8.9 Interpretation of red-infrared scattergram for 
pinyon-juniper landscape, aerial data. 
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Table 8.3 


Estimated Mean Subpixei Parameters 
Pinyon-Juniper Watershed, Aerial Data 


Parameter 

Estimated Value 

R m^RED^ 

15.4 volts 

W 

79.3 volts 

E ^ R gI^RED^ 

15.5 volts 

E[R gI (A, R )] 

41.3 volts 

VA R[R g ,(A RED )] 

21.4 volts 2 

v AR[R g ,(A IR )] 

64.3 volts 2 

•t 

30.25* 

1 

0.294 

mean fractional 
canopy cover 
of watershed 

0.23 

\ 

1300 meters 2 

Soil line equation 

yv = 172 



Similar reasoning allows one to neglect the shadow terms in (7.6). Further, 
since gg is small in this case, it follows from (7.9) that 

VAR[m] 2 VARjgj] ( 8 . 10 ) 

Substituting (8.10) into (7.6), neglecting the shadow terms, and rearranging yields 


VAR[R(A)] 2 E[gj] 2 + VAR[gj 


VAR[R gi (A)] + [R m (A) - E[R (A)]] 2 VAR[gj] 

( 8 . 11 ) 


Thus, the seven equations (6.17), (7.8), (8.10), (8.9), (8.11) (the latter two written 

for both bands) were solved simultaneously to obtain estimates of R f A ) 

m v RED' 1 

^m^lR^ E[gj], E[g g ], VAR[m] and VAR[gj] for the entire ensemble of pixels 
covering the watershed. Those results are provided in Table 8.3. 


Finally, in order to obtain estimates of fractional cover on a pixel-by-pixel 

basis, equations (6.7), (6.24), (6.19) (with g R (A) neglected) and (7.15) were 

®s 

combined to yield 


R(A ir ) - <*R(A red ) + [R m (A IR ) - R m (A RED )]m + y(l - mf +1 (8.12) 

Equation (8.12) was solved explicitly for m for each pixel, that is, for each paired 
observation (RCAj^p), R(A jr )) in the data set. The computer programs necessary 
for the above analysis are provided in Appendix D. 

The results of the analysis are graphically displayed in Figures 8.9, 8.10, and 
8.11 and Table 8.4. Figure 8.8 indicates that the theoretical canopy-soil 
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ESTIMATED SUBPIXEL CANOPY COVER 



ACTUAL SUBPIXEL CANOPY COVER 


Figure 8.10 Estimated subpixel canopy cover versus 
ground truth obtained from video 
for pinyon-juniper watershed. 
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SUBPIXEL ILLUMINATED SOIL 
OR SHADOWED SOIL 



SUBPIXEL CANOPY COVER 


Figure 8.11 Comparison of estimated values of illuminated 
and shadowed soil cover to the actual ground 
truth obtained from video for the 
pinyon-juniper watershed, aerial data. 


164 




Table 8.4 


Comparison of 

Estimated and Actual Fractional Canopy Cover 
for Pinyon- Juniper Watershed, Aerial Data 


Time of 


m 

Acquisition 

(hours) 

Act. 

Est. 

10.2924 

0.47 

0.41 

10.2942 

0.21 

0.29 

10.3078 

0.27 

0.20 

10.3113 

0.11 

0.17 


Si 

gs 

Act. Est. 

Act. Est. 


0.39 

0.50 

0.14 

0.09 

0.74 

0.65 

0.05 

0.07 

0.65 

0.75 

0.09 

0.05 

0.86 

0.79 

0.04 

0.04 
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reflectance model and the parameter estimates of Table 8.3 lead to a "triangular" 
interpretation of the actual red-infrared scattergram. That interpretation stems 
directly from the following pair of equations obtained from combining (6.24) and 
(8.7) and neglecting the shadow term, 

R ^REd) = m R rn^RED^ + ~ m )^ + R gj^RED^ (8.13) 

r < a ir> = m W + d - m )^‘ M A m ) (814) 

For instance, the apex of the triangle in Figure 8.8 occurs at full canopy cover 
(m = 1) in which (8.13) and (8.14) reduce to 

R ( A REd) = R m^RED^ ( 8,1S ) 

r (Air) = R- m (A IR ) (8.16) 


The base of the scattergram occurs for bare soil (m = 0) in which (8.13) and 
(8.14) become 

R ( A REd) = R gj^RED^ ( 8 - 17 ) 

R ( V = R g / V ( 8 - 18 ) 

and where (8.17) and (8.18) are related by the soil line equation given in 
Table 8.3. 

Lines of constant canopy cover indicated in Figure 8.9 are obtained by 

setting m constant in (8.13) and (8.14) and letting R (A) vary according to the 

% 

soil line equation. In a similar manner, lines of constant soil reflectance are 


166 



established by selecting one pair of Rg (\ ED ) and Rg (A IR ) (one point on the soil 

line), and letting m range from 0.0 to 1.0 in (8.13) and (8.14). Three such lines 
are drawn in Figure 8.9 representing three different soil background reflectances. 


Finally, Figures 8.10 and 8.11 and Table 8.4 show a comparison, for Five 
pixels, of the estimated fractional covers with the actual ground truth estimates 
obtained from the video. Figure 8.10 contains a plot of the estimated values 
canopy cover versus the ground truth values (s = 0.061). Figure 8.11 contains a 
graph of both the estimated theoretical curves of gj(m) and gg(m) and the ground 
truth values. The good agreement between estimated values and the ground truth 
in both Figures supports the applicability of this method for estimating subpixel 
fractional cover of semivegetated scenes, at least for the pinyon-juniper landscape. 


8.3 Pinvon-Juniper Watershed: Landsat Thematic Mapper Data 

This example investigates the fractional cover of the same watershed used in 
the previous case, except that Landsat TM data are used instead of aerial 
observations. Also, a different version of the inverse procedure is used due to the 
nature of the scattergram. 


8.3.1 Landsat Thematic Mapper Data 

The Landsat TM data used in this analysis were extracted from Scene ID 
Number Y504771733XO, obtained on June 21, 1985. The scene was purchased 
from EOSAT, Lanham, Maryland. The data were ordered with standard 
radiometric corrections that remove possible sensor error according to EOSAT 
procedures, and with an original pixel size of 30 meters. That size, the same as 
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tli6 raw data, was selected Id order to eliminate resampling, and hence unnecessary 
distortion to the DN values of the original data. 

The month of June was selected as it is normally a dry period of the year 
when there are few clouds and the atmosphere is clear. That dryness enhances the 
analysis due to the relatively high contrast between the reflectances of the conifer 
canopy and the ground, which contains bare soil and senesced grasses. However, 
since that dryness can also increase the aerosol count in the lower atmosphere and 
thus the diffuse radiance. It was decided to select an image in which those 
atmospheric effects were minimal. 

Since there was no practical means to estimate aerosol density over the 
Beaver Creek Basin for the archived TM images, the selection of the specific scene 
was made on the basis of an indirect and qualitative assessment of the cloudiness, 
image clarity, and soil moisture for the available TM scenes which were taken in 
the month of June. Cloudiness was assessed by comparing the microfiches of 
several scenes, made available by the EROS Data Center, Sioux Falls, South 
Dakota, and then selecting those which appeared to be the clearest. Next, the 
precipitation and climatic records of nearby meteorological stations for the months 
April through June were examined for each of those scenes in order to select the 
lowest precipitation, and thus by association, a low soil moisture for the period up 
to and including the time of acquisition. Although the above approach yielded 
perhaps the clearest of all available images for the month of June and provided 
some assurance that the soil was relatively dry, it did not provide any 
quantitative information on important properties such as optical thickness and 
diffuse radiance. Those quantities can best be obtained using ground-based 
instrumentation during the time of acquisition. 
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The Beaver Creek watershed was visually identified on the original scene and 
extracted using the General Image Processing Software (GIPS) developed by Peter 
Ford at the Center for Space Research, M.I.T. The watershed was located within 
a rectangular area covering 30 (vertical) by 29 (horizontal) pixels, or a total of 
870 pixels. 

8.3.2 Red-Infrared Scattergrams 

Two scattergrams from the original scene are plotted in Figures 8.12 and 
8.13, respectively, in terms of the satellite DN values. Figure 8.12 includes a large 
region covering 235 sq. km. containing several watersheds and a variety of 
vegetation types and densities. It possesses a typical triangular shape with a 
curved top and flat base, characteristic of semivegetated regions. The lower left 
portion of the scattergram, somewhat detached from the main part, can be shown 
using topographic maps to represent principally water bodies and regions with 
extensive shadows such as cliffs and gorges. The soil line was thus defined from 
the locus of pixels at the base of the major portion of the scattergram as shown in 
Figure 8.12. The mean and variance of the soil line are given in Table 8.5. 

The portion of the scattergram associated with only the small Beaver Creek 
watershed is outlined on Figure 8.12 and also plotted separately in Figure 8.13. 

That scattergram is located entirely within the lower portion of the large 
scattergram and does not possess a triangular shape nor a flat base. Thus, no soil 
line could be discerned on the basis of the pixels located within the Beaver Creek 
Watershed. 
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Figure 8.12 


Red— infrared scattergram of 235 sq. km region encompassing 
Beaver Creek Watershed, Landsat TM data. 
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Table 8.5 


Estimated Mean Subpixel Parameters 
of Beaver Creek Watershed, 
Landsat TM Data 

Solar Zenith Angle, 0 = 27.5 degrees 
Similarity Parameter, t| = 0.241 


Percent Diffuse Radiation: 

Q 

Quantity Units 


I'^red) dm 

0 

lUm) DN 

0 

m 

0.25 

Si 

0.70 

8s 

0.05 

R'm(^RED) DN 

40 

R’m &m) DN 

119 

E[R'gi(X,R£D)] DN 

69.5 

E[R gj(^iR)] dn 

56.8 

Scene Parameters 


E[DN(A.re D )] = 58.8 


E[DNa m )] = 69.5 


Soil Line Parameters 


E[DN(Xr£ D )] = 69.5 


E[DN(X m )] = 56.8 


Slope, a = 1.03 


Intercept, (3 = -15.4 



IQ 


2Q 

6 

12 

18 

6 

12 

18 

0.23 

0.21 

0.17 

0.72 

0.75 

0.79 

0.05 

0.04 

0.04 

30 

18 

1 

116 

114 

123 

63.5 

57.5 

51.5 

50.8 

44.8 

38.8 


VAR[DN(Xre D )] = 54.2 
VARtDNO.jR)] = 16.1 


VAR[DN(Xre D )] = 141.6 
VARtDNfXm)] = 161.3 
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8.3.3 Estimation of Subpixel Canonv Cover 

The inverse procedure used in this example was the same as that used for 
aerial case, except that the moment equations were written in terms of the DN 
values which include atmospheric effects. Further, the soil line, which could not 
be defined from the Beaver Creek scattergram, was assumed to be the same as 
that of the larger image. 

As in the aerial case, the Pinyon-Juniper landscape was assumed to consist 
of Poisson distributed spheres. The similarity parameter, rj, was calculated to be 
0.241 using Table 6.1, Item (iv), and an estimated solar zenith angle of 27.5 
degrees (after Iqbal, 1983, Chapter 1). The Sampling Scale Ratio for Poisson 
distributions, equation (6.23), yields 

Sp = 700/0.241A t >> 10 (8.19) 

and thus a large Sampling Scale Ratio was assumed, allowing one to use (6.24) in 
the analysis. 

The analysis assumed that the atmosphere was horizontally homogeneous, 
that the landscape was regionally homogeneous (no sharp contrasts in landscape 
reflectance), and that the vegetation reflectance was constant in both wavelengths. 
Further, since tj is small, the shadow terms were neglected as previously argued in 
(8.9) and (8.11). Under those assumptions, the expected value and variance of the 
entire ensemble of pixels in the scattergram can be written in terms of the 
satellite DN values using equation (6.31) and (6.32), respectively, or, 

E[DN(A,x)|m] = mR^(A) + gl E[R^(A,x)] + 1^(A) (8.20) 
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and 


VAR|DN(A,s)|m] = VA R (mir (A) + g,R' (A,*)' 


( 8 . 21 ) 


where R^(A) and R' (A,x) are "effective" reflectances of the vegetation and soil as 
®I 

previously used in Sections 4.3 and 6.3, and 1^(A) is the backscattered solar 
radiation term (or simply a calibration coefficient) characteristic of that scene. 

The soil line equations can be obtained by conditioning (8.20) and (8.21) 
along m = 0, or 

E[DN(A,x)|m = 0] = E[R' (A,x)J + U(A) (8.22) 

S I u 

and 

VAR[DN(A,x) |m = 0] = var[r'(A,x)1 (8.23) 

l Sj J 

It is noted that the calibration constant is included in the expected value of the 
soil line, and thus mean values of the soil reflectance can not be determined 
explicitly as in the previous examples. The variance equation can be expanded, 
analogous to (8.11), by inserting (8.10) into (8.21) and rearranging, or, 

VAR[DN(A)1 s [Elg,] 2 + V AR[gj]J VA R [R ' ( A)] + [RJA) - E[Rg (A)]j 2 VARfe,] 

(8.24) 

The number of unknowns include the fractional cover variables, m, gj, gg, 
VAR[gj], the reflectance quantities, R^(A red ), R^(A IR ), E[ R ' (A red )], 
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E[ R gj( A IR)]’ VAR i R g (^red^’ VAR [ R g (^ir)J* the atmospheric coefficients, 

>d< W “ d ‘d (Air), for a total of twelve. The equations available to solve the 
problem include i) the mean and variance of the observations, (8.20 and 8.24), 
written for both bands, ii) the mean and variance of the soil line equation, (8.22) 
and (8.23) written for both bands, and iii) the cover relationship gj(m) or (6.24), 
and (8.6) for a total of ten. Thus, the inclusion of the diffuse radiation terms 
results in two more unknowns than equations. 

There are several ways in which the diffuse terms can be estimated without 
actual atmospheric measurements. One method is to recognize that observations 
over areas of approximately zero reflectance (i.e. deep clear water bodies) consist 
principally of the diffuse radiance (Lillesand and Kiefer, 1987). A second approach 
is to use observations of two additional visible bands over pixels in which the 
reflectance is assumed independent of wavelength. By assuming an aerosol 
distribution, it may be possible to estimate optical depths and diffuse radiation by 
examining the relative intensity of the two bands (Liou, 1983). A third approach 
is to develop two additional independent equations using the moment analysis, 
such as the use of the cross-spectral covariance or two separate conditional lines. 

Since it was not possible to verify the estimated diffuse radiation using any 
approach, it was decided to simply assume a range of values and to compare the 
results as a function of those assumed values. In general, for optical depths of 
about 0.1 or less and surface reflectances of about 0.3, it can be argued that the 
backscattered diffuse radiation ranges from about 10 to 30 percent of the total 
radiation observed by the nadir-viewing satellite, depending on the wavelength. 
Scattering is likely to be greater at smaller wavelengths (i.e. blue) due to the 
combined effect of both Rayleigh and aerosol scattering. In the red and infrared 
regions, aerosol scattering is likely to contribute the most to the total optical 
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thickness (Liou, 1983). 

The above set of equations (8.20, 8.24, 8.22, 8.23, and 6.24) written, when 
applicable, for both wavelengths, were solved assuming that the backscattered 
solar radiation term, 1^(A), equaled 0, 10, 20, and 30 percent of the mean DN 
value of all the pixels in the watershed. For example, since the mean value in the 
red and infrared bands equaled 58.8 and 69.5 DN's, respectively, then for the 20 
percent case, the backscattered radiation was assumed to be 12 DN's. For 
simplicity, both bands were assumed to possess the same diffuse radiance, ignoring 
wavelength dependency. 

The solution to the above set of equations was found by minimizing the 
error between the theoretical moments and the actual moments in the same 
manner as for the aerial case. However, it was observed that in a few instances, 
two or three minimum values were obtained thus yielding two or three possible 
solutions. In all cases, the two additional possible solutions occurred at extreme 
values of m (very large, 0.90, or very small, 0.05) that were clearly inappropriate 
and thus they were not selected as the best estimate. That choice was supported 
by an a priori knowledge that the vegetation cover of the region was neither 
extremely dense nor sparse. However, future investigations would not necessarily 
benefit from such knowledge. 

The results of the analysis are provided in Table 8.5 and Figure 8.14. 

Table 8.5 indicates that the estimates of fractional cover are about 0.20, in good 
agreement with field observations (about 25 percent). The results are moderately 
sensitive to the relative magnitude of the diffuse radiation. The greatest value, 
m = 0.25, is estimated when no diffuse radiation is assumed. The estimate 
gradually decreases to m = 0.17 for the case with 30 percent diffuse radiation. 

Figure 8.14 shows the hypothetical data spaces of the 0 and 20% diffuse 
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Hypothetical data space 
20 % diffuse radiation 
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Figure 8.14 Hypothetical red-infrared data space for different amounts of 
backscattered diffuse radiation. 
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radiation cases, as compared to the actual scattergram. The hypothetical 3 pace 
for the 0% case was obtained in the same manner as the aerial case using (8.13) 
through (8.18), except that R(A) was replaced by DN(A), and the component 
reflectances were replaced by the effective reflectances. The data space for the 
20% case was obtained using the same procedure, except that the assumed values 
of diffuse radiation were added to those equations. The results indicate that the 
peak of the 20% case is located further from the soil line than the 0% case. Since 
the positions of the actual scattergram and soil line do not change, the estimate of 
fractional cover is directly related to the relative position of the hypothetical peak. 
As the peak moves further away from the soil line, the estimate of fractional cover 
decreases. 

An alternative way of plotting the 20% case is to use (8.13) through (8.18) 
as described above without adding the diffuse terms. The "adjusted" data space 
would be identical in size and shape as the 20% case plotted in Figure 8.14, but 
the entire space would be shifted 12 units downward and 12 units to the left. 

While the overall results of the analysis provide analytically reasonable 
quantities, they indicate the increasing difficulty in the ability to estimate 
fractional vegetation cover when the diffuse component becomes large. Of 
particular concern is the rapid decrease in the estimated vegetation reflectance 
with increasing diffuse radiation. When the diffuse component is at 30 percent, 
the vegetation signal is only (1)(0.17)/18, or about one one-hundredth of the 
magnitude of the diffuse term. The signal to noise ratio is thus very small and 
even relatively minor perturbations in atmospheric effects would violate the 
assumption of horizontal homogeneity. Thus, accurate inversion for vegetation 
properties under conditions of moderate diffuse radiation seems unlikely. 
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Chapter 9 


SUMMARY 

9.1 Principal Conclusions 

The research in this report has demonstrated that it is theoretically feasible 
to estimate spatially-variable bulk properties of semivegetated landscapes at 
subpixel scales for optically-thin atmospheres using only one set of multispectral 
observations without ground truth, at least for a limited range of landscapes. The 
approach relies on the physically-based conceptualization of landscapes as 
stochastic-geometric reflecting surfaces, which can possess variability in both the 
geometry of the shape and spatial distribution of the plants, as well the vegetation 
and soil background reflectance. The degree to which subpixel parameters can be 
retrieved depends on several factors including knowledge of the structure of the 
landscape, the number of landscape variables, the magnitude of the Sampling Scale 
Ratio, and the ability to identify groups of pixels within the red-infrared 
scattergram which possess common attributes. 

An important feature of the inverse procedure is that it takes advantage of 
the multispectral nature of the data by solving equations associated with both the 
red and infrared wavelengths simultaneously. It thus extends the work of others 
(Otterman, 1984; Li and Strahler, 1985) who have inverted geometric models using 
only one band and assumed reflectances. The methodology offers a 
physically-based alternative to current practices which are highly empirical. 

The reflectance model and inverse technique are primarily applicable to 
regional scale hydrologic investigations where the parameterization of numerous 
plant and soil properties is not feasible nor of practical importance at such large 
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scales. By absorbing the variability of such properties into a few bulk plant and 
soil variables, an inherent tradeoff is made between the amount of physical detail 
which can be modeled or estimated, and the size of the region which can be 
investigated. The technique relies on the existence of a large number of pixels 
possessing a wide range of soil and vegetation. Thus, when only a few pixels are 
available, or when the entire scene is homogeneous, a different version than those 
presented above should be considered. 

The inverse method has been tested only on idealized simulated scenes and 
one conifer watershed using both aerial and satellite data. Good results were 
achieved in the case of Beaver Creek despite some major assumptions including 
the Poisson distribution of the trees, constant vegetation reflectance, especially for 
the infrared band, similar soil reflectance for the soil line and the semivegetated 
areas, and the neglect of shadow contributions for small 77. Further, the procedure 
had to be adapted to each case based on the shape of the scattergram and the 
limited knowledge of the landsurface. While those algorithms worked well for 
both the idealized and actual cases, their general applicability to other 
semivegetated landscapes is unknown. Thus, further testing is warranted on a 
wide range of other types of semivegetated landscapes in order to validate and 
improve upon the methodology presented in this report. 

The landscape reflectance simulation model developed in this research has 
been shown to be an effective mechanism for investigating the sensitivity of 
landsurface variability on the behavior of multispectral data acquired at scales 
representative of current satellite pixels. That feature provides a useful alternative 
to the toilsome and expensive task of understanding variability in actual scenes by 
obtaining simultaneous ground truth for a large number of pixels. By sequentially 
altering different variables into the simulations, valuable insight on the 
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multispectral behavior of actual images can be obtained. Modeling of the 
landscape is facilitated through the introduction of the non-dimensional similarity 
parameter, 77, which generalizes the results without constraining them to any one 
geometric shape or solar angle. 

A principal finding from the simulations has been the recognition of the 
correlation which develops among fractional canopy cover, shadow, and illuminated 
soil background with increasing pixel scale. The analytical formulation of the 
gj(m) relationships for different spatial distributions (Figures 7.8, 8.4, and 8.11), 
based on either deterministic or statistical reasoning, has facilitated the solution of 
the inverse problem by eliminating one or two unknown parameters. That 
correlation has been shown to be a principal mechanism that contributes to the 
evolution of the tasseled cap of red-infrared scattergrams of semivegetated 
landscapes. 

The moments of the reflectance equations have been expanded in terms of 
the moments of the individual variates of the stochastic-geometric reflectance 
model for the purpose of solving the inverse problem: The estimation of subpixel 
parameters given only the red— infrared scattergram and limited assumptions on the 
structure of the scene. The inverse procedure involves equating those analytical 
moments to the actual moments of the image, and solving the equations 
simultaneously, without the need for ground truth. 

Knowledge of the relationship between the physical structure of the 
landscape and the shape and structure of the scattergram has been shown to 
facilitate the inverse problem in at least two manners. First, it provides a 
mechanism for identifying pixels with common attributes, especially for cases with 
large Sampling Scale Ratios. Second, it allows the formulation of additional 
moment equations, conditioned on those common attributes, which are often much 
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simpler than the general moment equations required for cases with small Sampling 
Scale Ratios. 

The exact formulation of the inverse procedure depends on several factors 
related to the structure of the landscape and to the interpretation of the 
red-infrared scattergram. There is no single recipe that can be listed that 
accommodates all situations. However, several general steps common to the 
solution of most of the inverse examples given in this report are as follows: 

i) All the radiometric observations for a particular region are plotted in the 
red-infrared data space. 

ii) A narrow band of pixels lying at the base of the scattergram is selected 
as representative of the soil background line, and the reflectance moments of that 
ensemble of pixels are computed. 

iii) The solar zenith angle is computed based on the time of overpass. 

iv) If possible, an assumption is made on the bulk geometric shape of the 
plants on the landscape. The similarity parameter, tj, is then calculated based on 
the plant shape and solar zenith angle. If the shape is unknown, then r/ may have 
to be included as an unknown in the analysis. 

v) The Sampling Scale Ratio is determined based on the scale of the pixel, 
the similarity parameter, and an assumed value of the horizontal scale of the tree. 
This calculation is only an order of magnitude estimate and an exact value of the 
tree size is not required. 

vi) For small Sampling Scale Ratios, the inverse problem is solved using the 
full set of moment equations. Their number and complexity will depend on the 
number of assumptions one is willing to make on the structure of the landscape. 

In some examples, the moment equations can be simplified by neglecting relatively 
small terms, such as the shadow terms when rj is small. In general, only the 
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moments of the fractional cover types for the entire ensemble of pixels can be 
retrieved for small Sampling Scale Ratios. 

vii) For large Sampling Scale Ratios, several conditional moments can be 
written for pixels with common attributes, identified through ones knowledge of 
the shape and structure of the scattergram. In general, for cases in which 
vegetation reflectance can be assumed constant, pixels of equal vegetation amount 
will orient themselves parallel to the soil line. In such cases fractional cover can 
be estimated on a pixel— by— pixel basis. 

The case studies indicate that the idealized tasseled-cap scattergram is not 
always realized due to limited combinations of soil and vegetations properties. 
However, inversion can still be achieved if a representative soil line is identified. 

Finally, the analysis of the Thematic Mapper data indicates that subpixel 
canopy cover can be estimated on a pixel-by-pixel basis without specifying the 
absolute magnitudes of the soil and vegetation reflectances, even when there exists 
a small backscattered diffuse component. That conclusion is limited to cases with 
horizontally-homogeneous atmospheres and regions of low contrast. Preliminary 
results indicate that neglecting the diffuse radiation component tends to 
overestimate the mean fractional cover of the region. 

9.2 Future Research 

There are several directions in which future research can proceed. First, 
there is a need to understand the sensitivity of the present conclusions to the 
various assumptions required for inversion. That can best be achieved with the 
aid of the simulation model. The model can be extended to include other factors 
including different spatial distributions of plant spacing and soil reflectance, 
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atmospheric effects, and other similarity values. Inversion techniques can then be 
applied to the new set of scenes. A further extension which would improve the 
flexibility of the model is the incorporation of topographic effects, for instance, in 
conjunction with U.S.G.S. digital elevation maps. 

A more rigorous analysis of the effects of diffuse radiation, both surface 
reflected and solar backscattered, is warranted in order to understand its effect on 
both the structure of the scattergram and the estimation of subpixel properties. 

For the inverse procedure, it might be possible to include the diffuse terms as 
unknowns into the analysis. The limits to which the coupled 
landsurface-atmosphere radiation model is applicable needs to be defined in terms 
absolute values of the principal surface and atmospheric properties, such as optical 
depth, landsurface reflectance, and solar zenith angle. 

An important application of the reflectance model is the understanding of 
the physical basis of common vegetation indices. Although the vegetation indices 
are very empirical, they are nonetheless widely used by scientists for the 
assessment of vegetation amount. The simulation model can be used to generate 
common indices to investigate the sensitivity of subpixel variability on the shape 
of vegetation indices, in a manner similar to that used in the present report for 
the understanding of red-infrared scattergrams. 

The original motivation for the development of the present research was to 
define landscape properties necessary for the application of the equilibrium 
hypotheses noted in Section 1.1, specifically with regard to estimating soil 
hydraulic properties. The analysis of the Beaver Creek vegetation in this report, 
together with additional data already collected (Jasinski, 1987), provide the 
necessary input to complete that study. 
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Appendix A 


Red-Infrared Data Spaces 
of Simple Hypothetical Semivegetated Scenes 


185 


</T 

<L> 


i-^s rr\ 



in 


c 

C 

ro C 

CQ i— i 


186 


Band 2 Intensity 

Figure A.l Hypothetical red-infrared space of constant vegetation reflectance landscape 


Effect of increasing R4 
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Effect of 
Increasing R4 



188 


Band 2 Intensity 

Figure A. 3 Effect of increasing soil reflectance on hypothetical red-infrared space 


Original Data Space 
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Band 2 Intensity 

Figure A.4 Effect of shadows cast by cones on hypothetical red-infrared space. 








Original Data Space 
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Band 2 Intensity 

Figure A. 5 Effect of variable vegetation reflectance on hypothetical red-infrared space. 


Appendix B 


Analysis of La ndsat 2 Multi spectral Scanner Data for the Tans gtudx Area 

During the early stages of this research, several linear regressions were 
conducted between formulas using Landsat 2 MSS data and fractional cover 
estimated from aerial photography, for a region centered over Taos, New Mexico. 
That analysis was conducted prior to the development of the canopy-soil 
reflectance model and inverse procedures presented in the main portion of this 
report. The results of those linear regressions are presented in this appendix. 
Details of the work were reported by Jasinski and Eagleson (1986). 

It is noted at the outset that the correlation coefficients of the regression 
analyses were low and the results of the regressions were considered inconclusive 
due to a variety of reasons. The primary reasons included i) uncertainty in the 
quality of the MSS data, which had gone through several preprocessings including 
at least two resamplings, ii) problems in registering a given Landsat pixel to a 
particular location on the aerial photograph, and iii) difficulties in estimating 
fractional cover using the color aerial photographs. Nonetheless, some insights 
were gained and a summary of the regression analyses is provided below. 

B.l Site Description 

The Taos Study Area is outlined in Figure B.l. The land includes a wide 
variation in surface relief, ranging from flat plains to rolling foothills, to detached 
high ridges. Elevation ranges from 6,000 to 10,000 feet. Vegetation tends to 
follow the topography. The lower flats are covered with blue grama and 
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Figure B.l Location Map, Taos Study Area. 
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wheatgrass grasslands, and snakeweed, rabbitbrush and sagebrush shrublands. 
Pinyon-juniper woodlands are found in the rolling foothills. At the higher 
elevations, there is ponderosa pine, spruce, fir and aspen. Percent cover ranges 
from nearly 0 to 100 percent, with the majority of the area 40 to 60 percent 
covered. At least two trends in percent vegetation cover can be readily observed. 
They are, first, a decreased percent vegetation cover with decreasing altitude, and 
second, a less dense cover on south-facing slopes compared to north facing slopes 
at the same elevation. 

B.2 Aerial ^nd Satellite Data 

The database consisted of Landsat MSS data, and 1:3000 aerial photographs, 
supplied by the Bureau of Land Management, Branch of Remote Sensing, Denver, 
Colorado and the Bureau of Land Management, Taos Resource Area Office, Taos, 
New Mexico (Work, 1983). 

Lgjnjsat Data . The Landsat scene used for this analysis was derived from an 
original Landsat MSS scene, Number 21608-16562, on June 18, 1979. The scene 
included some preprocessing by BLM in addition to that routinely supplied by the 
EROS Data Center on original CCTs. The processing consisted of 1) the removal 
of certain radiometric and electronic anomalies known as line drops and banding 
by filtering, 2) the removal of minor geometric distortions which were inherent in 
the data, and 3), the registration of the Landsat data to a Universal Transverse 
Mercator (UTM) map projection using a resampled 100 meter square pixel. Since 
the regression procedure worked on a pixel-by-pixel basis, correct registration was 
of paramount importance. Landsat data were fitted to the UTM grid by visual 
inspection through the use of color slides of Landsat segments projected directly 
onto USGS topographic maps. 


193 



Aerial Photography . Approximately eighty color aerial photographs at about 
1:3000 nominal scale were borrowed from the Bureau of Land Management, Taos 
Resource Area Office, Taos, New Mexico for the current study. Those were taken 
on June 16, 1981 using a relatively low flying aircraft with a nine inch square 
format and a six inch focal lens. 

Photographs were selected to represent a broad range of vegetation cover and 
to exclude agricultural and urban areas. Because of the random nature of the 
photograph locations, the eighty photographs were distributed over twenty 
different USGS quadrangles. 

Photograph analysis included several steps. First, photographs were visually 
registered to the UTM grid by comparing topographic features of the photograph 
to those of the USGS map. Next, the photograph was divided into pixels 100 
meters square using a clear overlay and the center eight pixels were selected from 
each photograph. Each pixel at 1:3000 scale was about 1 to 1-1/2 inches square 
and contained a random vegetation cover interspersed with soil background. 

Fractional vegetation cover for each pixel were analyzed using an image 
analyzer connected to a video camera. Percent cover was determined by selecting 
for each pixel the threshold "grey level" associated with only the vegetation cover 
and then computing the total area below (darker than) the threshold level. The 
procedure worked satisfactorily for pixels which contain distinct vegetation and soil 
characteristics. The error of the canopy cover estimate for such cases, which 
represent roughly one half of the over 100 pixels analyzed to date, was several 
percent. For pixels, containing dark soils or significant shadows, error was 
estimated to be roughly ±10 percent. Roughly twenty percent of the pixels 
analyzed fit into the latter category. 



B.3 Regression with Normalized Difference Vegetation Indev 

This analysis consisted of regressing the NDVI with ground truth obtained 
from the aerial photographs. Two variations of this approach were tested. A 
total of 116 pixels were used. The first variation involved using the NDVI defined 
in terms of actual integer DN values instead of reflectances. The second approach 
used actual radiances computed using conversion factors described by Markham 
and Barker (1986). The results are provided in Table B.l below and shown on 
Figures B.2 and B.3. 


Table B.l 

Normalized Vegetation Index 
versus Percent Cover 

NVDI Variation m 


B2 


2*DN 4 - DN 2 

VI DN = 5 ? DN 4 + DN 2 m = 199 VI DN + 095 °- 61 

VI r = K^- + - 100 m = 3.06 VI R - 111 0.58 

The results indicate that for both variations, about 60% of the change in 

NDVI can be explained in terms of percent vegetation cover. 

B-4 Regression Using Direct Beam Emiatinn 

Assuming that the landscape consists of only two cover types, soil and 
vegetation, equation (4.10) can be rewritten, 
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Figure B.2 







»y cover. 




m 


= L' - R g]/[ R m - R g 


(B.l) 


For constant reflectances, m is thus linearly related to the observed direct beam 
radiance or the DN values. With this in mind several linear regressions were 
carried out with m as the dependent variable and the DNs as independent 
variables. The results are shown in Figures B.4 and B.5 and summarized in Table 
B.2, 


Table B.2 

Summary of Linear Regressions 
Direct Beam Equation 


Recession Eauations 

B2 

m = 126 - 199DN2 

0.53 

m = 377 - 138DN4 

0.10 

m = 120 - 197DN2 cos 0 

0.53 

m = 143 - 189DN2/cos0 

0.42 


where DN2 and DN4 represent MSS bands 2 and 4, respectively. 

Variations of the above direct beam equation include accounting for changes 
in zenith angle due to topographic slope. For surfaces with average slope /?, the 
equation (B.l) can be rewritten 


m 


^ - RJ/[R m - R 

L g J/l m ' 


gj 


(B.2) 


Zenith angle effects can be included by 
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Figure B.4 
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m 


(B.3) 


Lcos 3 
L ' cos 0 



where 0 is the ground slope, measured from USGS 7.5-minute topographic maps 
Zenith angle was computed using the following formula from Iqbal (1983), 


Cos 6 = (sin^ cos /? - cos <j> sin/? COS7) sin£ 

+ (cos^ cos /? + sin 0 sin/? COS7) cosS cosu 
+ (cos<? sin/? sin7 sinw) 

where 

6 = solar declination 
u ) = hour angle 

7 = surface azimuth angle 
/? = average slope of pixel 


Solar declination and hour angle were estimated from the time of Landsat 
overpass. 

The results of the linear regression including the ground slope and zenith 
angle corrections were also poor, as indicated on Figures B.6 and B.7, and on 
Table B.2. Although part of the explanation for the poor correlation may simply 
be due to the bidirectional reflection characteristics of the soil, a more likely 
explanation may also simply be the inaccuracies introduced by measuring small 
distances off the topographic maps. At the 1:24000 scale, pixels are less than 
0.2 cm 2 in area and small inaccuracies in measurement or pixel registration can 
cause serious error in the regression analysis. 
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B.5 Regression with Kauth-Thomas Indices 

The Kauth-Thomas greenness and brightness indices were computed and 
then regressed with actual percent cover obtained from the aerial photographs. 
The results are shown in Table B.3. They indicate, contrary to expectation, that 
brightness appears to explain more of the variation in m than greenness. 

Table B.3 

Regressions with Kauth-Thomas Indices 


Index 

BI 

Greenness, GU 

0.24 

Brightness, BI 

0.39 

regression analyses yielded, 

m = 10.47 + 0.14 GI 

(B.4) 

m = 82.82 - 0.48 BI 

(B.5) 


B.6 Multiple Linear Regression 

Multiple linear regressions were carried out with the same data set as in 
previous cases with percent cover as the dependent variable and the MSS band 
observations as independent variables. The two cases examined were m vs. DN1 and 
DN2, and m vs. DN1, DN2, DN3 and DN4. Once regression coefficients were 
obtained, theoretical percent cover obtained from the multiple linear regression 
analysis was regressed with actual percent cover in order to compare correlation 
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coefficients with other methods. The results for the second case are shown on Figure 
B.8 (for regression with four bands) and summarized on Table B.4. 


Table BA 

Results of Multiple Linear Regression 

Regression Equation j ^2 

m = -2.25 DN2 + 0.70 DN4 + 74.97 0.53 

m = -2.07 DN1 - 0.62 DN2 + 0.20 DN3 + 0.63 DN4 + 72.25 0.58 

As expected, there is negative correlation with the visible bands and positive 
with the near infrared. It is also noted that the addition of bands 1 and 3 only 
contributes an increase of 0.05 in R 2 . 

B - 7 Regression Using Linear Distances in the Red-Infrared gcattergram 

The fractional cover of a given pixel estimated from the aerial photographs, 
m g , was regressed with estimate of m based on the pixel's location in the 
red— infrared scattergram plotted in Figure B.9 (An expanded version of Figure 
1.6). The procedure consisted of the following: 

1) All data points within the segment were plotted. 

2) Envelope lines are drawn along the three sides of the triangular data 
space. The soil line was drawn as a straight line emanating from the origin. (For 
the assumption of no shadows and constant vegetation reflectivities, all sides of 
the triangle must be drawn straight.) 

3) Along the soil line, m was assumed equal to zero. Likewise, at the top 
of the triangle, m was assumed equal to one. 
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Actual Percent Green Canopy Cover, m 

Figure B.8 Percent cover from multiple linear regression, m , versus actual percent cover, 




Band 2 Diqital Number 

Figure B.9 Data plot of MSS Band 4 versus Band 2. 





4) For constant vegetation reflectivity, m was assumed linearly related to 
the distance between the top and base of the triangle. For example, for a point 
exactly halfway between the top and the base, m was assumed equal to 50% 
cover. 

The graphical results are shown in Figure B.10, which for clarity, includes 
only the data points with ground truth, and not the full scattergram. The R 2 
value resulting from the regression of m with m was 0.34. 

O 
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Appendix C 


Derivation of Equations 

This appendix describes the derivation and assumptions of the following 
equations presented in the main text: 

Cl.) General Reflectance Moment Equations 
Cl.l) Mean Reflectance 
Cl. 2) Variance of Reflectance 
Cl. 3) Cross-Spectral Covariance 
Cl. 4) Spatial Covariance of Reflectance 
C2.) Reflectance Moment Equations: Case V, Method 2 
C2.1) Mean Reflectance 
C2.2) Variance of Reflectance 

C3.) Geometric Similarity Formula for Poisson Distributions 
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Cl. General Reflectance Moment Equations 


The reflectance equation is 

R(A,x) = l fj(x) Rj(A,x) (4.1) 

i 

where the terms are defined in Section 4.1. 


Cl.l Mean Reflectance 


Considering all the terms in (4.1) as random variables, the mean reflectance 
in its most general form is 

E[R(A,x)] = £ E[fj(x) Rj(A,x)) (Cl.l) 

i 

Using the identity 

E[xy] = COV[x,y] + E[x] E[y] (Cl. 2) 

equation Cl.l can be expanded 


E[R(A,x)] = l E[fj(x)] E[Rj(A,x)] + COV^x), R i (A,x)]| (C1.3) 


Assuming that the fractional cover, fi(x) is independent of the reflectance of 
that fractional cover, Rj(A,x), the above equation becomes, 


E|R(A,j£)] = l E[fj(jt)] E[R i (A, i )] (C1.4) 

i 

which is the same as (6.1). 


Cl. 2 Variance of Reflectance 

The variance of (4.1) in its most general form is 
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var[R(A,x)] = l VAR[fj(x) R i (A,x)] 


+ II cov [l f jCs) Hj(A,x)], [fj(x) Rj(A,s)] 
i#j " J 


(C1.5) 


If fj(x) and Rj(A,x) are independent, the first term on the right hand side of the 
above equation can be written 


I VAR[fix) R i (A,x)]= ^[E[fj(x)] 2 var[Rj(A,x)] 
i i 

+ E[Rj( A, x)] 2 VAR[fj(x)] 

+ VAR [ f j(-)l VARfR^A^)]] (Cl. 6) 

If it is further assumed that the reflectances are independent of each other, the 
second term on the right hand side of (Cl. 5) becomes 


II cov 

Mj 


[ f j(x) Rj(A,x)], [f.(x) Rj(A,x)] 


= 11 R j( A ’*) R j( A >*) COV[fj(x), fj(x)] 

i#j 

Combining Cl. 6 and Cl. 7 yields the variance equation given in (6.2). 


(Cl- 7) 


Cl. 3 Cross-Snectral Covariance 

The cross-spectral covariance can be expanded using Cl. 2 as 
COV[R(A,x), R(A 2 ,x)] 

= E[R(A r x) R(A 2 ,x)] - E[R(A r x)]E[R(A 2 ,x)] (Cl.8) 
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Inserting (4.1) into (Cl. 8) yields 


COV[R(A r x), R(A 2 ,x)j 


= E 


[ l fj(x) Rj(A,,s)] [ l fj(x) Rj(A 2 ,x) 


- e[ l fj(x) RjfAjjt)] e[ l {j(i) R.(A,i)' 


(Cl. 9) 


Assuming that the fractional covers, fj(x), are independent of the reflectances, 
Rj(A,x) yields 

COV[R(A 1 x), R^)] 

= n E l R i( A r i) Rj(A 2 .s)] E[f,U) fj(x)] 
i j 


- II E [ R i< A r i) EI R j(A 2 ,x] E[f.(i)l E[f.(*)] 


(C1.10) 


i J 


Applying (Cl.2) to E[fj(x) fj(x)] and to E[Rj(A 1 ,x) Rj(A 2 ,x)] and inserting into 
(C1.10) yield 


COV[R(A 1 ,x), R(A 2 ,x)] 

= n [ c °v[ R j( A r x). Rj(AjS)] + E[R.(A lJt )] EIRjCAjjs)]" 
[cov(f.(a),f j ( 2[ )] + E[f,(*)] ElfjU)]] 


1 J 


- n E i R i( A i^)j E i R j( A 2’-)i E [ f i(i)) E i f j(s>] 


(Cl. 11) 


1 J 


Cancelling terms of opposite sign yields 
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COV(R( A R(A 2 ,x)) 

= II [cOv|R j (A 1 ,x), Rj(A 2 ,2c)] COV^x), 
i j 

+ COVlRjfAjji), R j (A 2 , 2£ )| Effjti)) E[fA(x)J 



+ E [ R jU r x)] e [Rj(A 2 ,x)J COV(fj(x), fj(x)] 

(Cl. 12) 

When i = j, then 



and, 

fj(x) fj(x) = fj(x) 

(Cl. 13) 


COV[fj(x), fj(x)] = VAR[fj(x)] 

(Cl. 14) 

Assuming there is no cross-spectral covariance between the reflectances of different 
cover types, then 


COV[R i (A 1 ,x), Rj(A 2 ,x)] = 0 for i t j 

(Cl. 15) 

Separating the summation in Cl. 12 into portions representing i 
j, and inserting Cl. 13, Cl. 14, and Cl. 15 yields, 

= j and i i 


COVjRfAj.x) R(A 2 ,x)] 

= I [E[fj(x)] 2 COv[Rj(Aj,3j), R i (A 2 , J! )] 
i 

+ VAR[fj(x)] ElR^x)} E^A^x)]] 

+ I I E [ R i( A !^)] E [Rj(A 2 ,x)] COV[fj(x), fj(x)] (Cl. 16) 

i^j 

which is the same as that given in (6.3). 
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C1.4 


The spatial covariance of the total reflectance between two pixels located at 
x and x' can be written 

COV x [R(A,x), R(A,x')] = COV x [ l fj(x) Rj(A,x)j , fj(x')Rj(A,x')l (Cl. 17) 

1 j 

where i and j represent cover types. The above can be expanded 

COV x (R(A,x),(A,x-)I = E [ l fjU) Rj(A, 2£)] [ l fjfa' ) RjfA.a')] 

i j 

- e[ l fjte) Rj((Aj£)] e[ l fj(r) Rj((A,x-)] (Cl. 18) 
> j 

The first term on the right hand side of (Cl. 18) can be rewritten 

E [ l fjtx) Rj(Ajt)] [ J fj(i') Rj(Ajt')] 
i j -I 

= E [ 1 1 *j(x) f j(2c') Rj(A,x) Rj(A,x')j (Cl. 19) 

i j 

Assuming the fractional covers are independent of the reflectances, then 
(Cl. 19) becomes 

e[ n f i(*) f / r) R j< A J£') 

» j 

= 1I E [ f i(-^ E [ R j(A^) Rj(A,x') (Cl. 20) 

i j 

The same assumption allows the second term on the right hand side of 
(Cl. 18) to be written 


215 



(Cl. 21) 



= l E[fj«] E[Rj(Ajj)] l E[f (x')] E[R (A,r)] 

i j 1 

Equations (C1.20) and (Cl. 21) can be rearranged and combined to yield 
COV x [R(A,x), R(A,x')] = 

1 1 [ElfjW fj(x')] E[Rj(A,s) R (A,x-)] 

i j 

~ E[fj(x)] E[fj(x')j E[Rj(A,x)] E[R j (A,x')]] 

Using the identity given in (C1.2), 

E[fj(x) fj(x')) = COV x [f i (x), fjtjt')) + Effjfjc)] E[f.( S ')J 

E[Rj(A,x) Rj(A,x')] = COV x [Rj(A,x), R.(A,x')j 

+ E[Rj(A,2>] E[Rj(A,x')] 

Inserting (CI.23) and (C1.24) into (C1.22) yields 
C°V x [R(A,i), R(A,x')] = 

l l [{COV x [fi(x). f.(r)] + E[f.( S ) f ( S -)1) • 
i j 

{COV x [R.(A,x), Rj(A,x')] + E[Rj(A,x)] E[R.(A,x')]} 
- EffjCx)] E[f.(x')] E[R.(A,x)j E[R.(A,x')] 
which is the same as (C6.4). 


(Cl. 23) 


(Cl. 24) 


(Cl. 25) 
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C2. Reflectance Moment Equations for Case V - Method 2 


The reflectance moment equations for Case V, Method 2, can be derived 
from equations (6.1) through (6.4). The notation is changed to that of equation 
(5.1), or 


R(A,x) = m I R m ^(A,2t) + m g R m (A,x) 

+ Sl R gi ( A ’*) + gg R gg ( A >s) 

where i = 1 (in equation (4.1)) designates illuminated vegetation, i = 2 designates 
shadowed vegetation, i = 3 designates illuminated soil background, and i = 4 
designates soil background shadowed by vegetation. 


C2.1 Mean Reflectance 


In this example m g = 0, and R (A) and R (A) are assumed constant 

1 ®s 


throughout space. The mean reflectance equation is thus 


E[R(A,x)j = E[m] KJX) 4- E[gJ R gg (A) + E[g,] E[R (Ajt)] (C2.1) 

Cl.2.2 Variance of Reflectance 

In this case it is recognized that the variance of all the constant reflectance 
terms equals zero, or 

VAR[R m (A)] = VAR[R gg (A)] = 0 (C2.2) 

The variance equation in (6.2) becomes 
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var[R(A,s) 1 = E[g ] 2 var[R (A,x)J 

6 I 

+ R m ( A ) 2 VAR[m] + E[R g ^(A,x)] 2 VARfgj] 

+ VAR[g ] VAR[R (A, 20] 

h 

+ R (A) 2 VAR[g ] 
g s s 

+ 2E f R m^ A ^ E [ R gi ( A ’S)] covfm^] 

+ 2E[R m (A)] E(R (A)] COV[m,g s ] 

+ 2E [ R g ( A ^)] E ( R g ( A )] COV[g r g s ] (C2.3) 

I s 

In order to reduce the number of unknowns, all the covariance terms will 
be expressed in terms of COv[m,g ]. For instance, 

COV[m,g I ] = COV[m, 1 - m - g g ] 

= E[m(l - m - g g )] - E[m] E[1 - m - g g ] 

= E[m - m 2 - mg ] - E[m] + E[m] 2 + E[m] Efgl 

= E[m] - E[m 2 ] - E[mg ] - E[m] + E[m] 2 + E[m] E[gl 

= -[E[m 2 ] + E[m] 2 J - [E[mg g ] - E[m] E[g g ]J 

= - VAR[m] - COV[m,g s ] (C2.4) 

Following a similar expansion as above for COv[g r g s ], and inserting those 
covariance expressions into (C2.3) yields the variance expression given in (7.6). 

The cross— spectral covariance, given in (7.7), can be obtained in a similar 
manner as for the variance. 
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This section derives the geometric similarity equation for Poisson 
distributions given in (6.17). Assume that a flat surface of area A p is covered 
with a random number of flat two-dimensional figures of arbitrary shape and size 
with mean area a. Assume that the figures can overlap one another, and that the 
centers of the figures are Poisson distributed in space with density p. The 
expected area that is covered, A c , can be shown to be (Kellerer, 1983), 

E[A c |p] = A p [l - exp(-a^)] (C3.1) 

letting f = A c /Ap then 

E[f|/>] = [1 - exp(-ap)] (C3.2) 

for trees or plants with projected mean area A t , the expected fractional area 
covered is 

E[m|/>] = [1 - exp(-A t />)] (C3.3) 

when the mean shadow cast by an individual tree, S t , is defined, 

S t = 7?A t (C3.4) 

then the expected fraction of plant and shadow is 

E[( m + g s ) I p\ = [i - exp[-/>(A t + S t )J 

= [l - exp [-p(r) + 1) A t j (C3.5) 

Assuming that the plant canopy overlaps the shadow (that is, there is no 
shadowed canopy), then 

E[m s ] = 0 
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From (5.7), 


E[g,J = 1 - [E[m] + E(g g ] 
Inserting C3.5 into the above, 

E[gj] = exp[~ P(*l + 1)AJ 
= [expf-Mj]^ 1 


Inserting C3.3 into the above yields 



which is the same as equation (6.17). 


(C3.6a) 


(C3.6b) 


(C3.7) 
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Appendix D 


Computer Programs 


Dl.) Canopy-Soil Reflectance Simulation Model 
D2.) Inverse Procedure for Case II 
D3.) Inverse Procedure for Case V, Method 1 

D4.) Inverse Procedure for Case V, Method 2: 

i) Estimate Soil Background Cover 

ii) Estimate rj 

D5.) Inverse Procedure for Beaver Creek: 

i) Estimate Bulk Parameters for Entire Scene 

ii) Estimate Bulk Parameters for Each Pixel 
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CANOPY-SOIL REFLECTANCE SIMULATION MODEL 


C 

C 

C 

C 

C 

c 

c 

c 


BY 

MICHAEL F . JASINSKI 

PARSONS LABORATORY FOR WATER RESOURCES AND HYDRODYNAMICS 
M.I.T, Room 48-212 
Cambridge, MA 02139 
Tel. 617-253-5483 


c 

consisting < 

c 

reflectance 

c 

The output : 

c 

(plotter) . 

c 

Version 14: 

c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 


c 

Variables : 

c 

XI, X2 - 

c 

RNODE - 

c 

RAGNODE- 

c 

PCTV 

c 


c 

PCTVAG - 

c 

PCTS 

c 

PCTSAG - 

c 

THGT 

c 

TREF 

c 

ETREF - 

c 

ETHGT - 

c 

TDIA 

c 

VARREF - 

c 

NVI 

c 

MREF 

c 

IBAND - 

c 

IAG 

c 

ISHAD - 

c 


c 

IELEV - 

c 


c 

IVREF - 

c 


c 

IVDIA - 


sparse 


canopy 
l. Total 


files are "CAN14 .OUT" (printer) and "C14BPL.DAT" 

computes visible and infrared reflectances, 
percent canopy, and percent shadow for four 
levels of aggregation (1,5,10,30). Includes 
calculation of percent canopy for 
each pixel. Includes soil variability determined 
from Turning Bands Model (calculated previously) 

For soil reflectance (ASSIGN SOILREF.OUT = FOR002) . 
Assumes relationship between tree vis/ir reflectances 
Increases dimensions to 150x150. 

Includes inhibitory field, ie.e* nonoverlapping 
°f trees . Trees positioned at center of one meter 
pixels. Includes shadows as an option. Version 14 
reorganizes the sequence of calculations of 
reflectance and shadows. If elevations are to 
be read then ASSIGN ELEV0_05.OUT FOR001. For 
shadows to be read, ASSIGN VEGREF . OUT FOR003. 

tree center coordinates 

actual reflectance at a given node, finest 
aggregated reflectance at a given node 
actual percent vegetation at a given pixel 
resolution (range 0-1.0) 

aggregated percent vegetation at a given node 
percent shadow in a given pixel (range 0 - 1.0) 
aggregated percent shadow 
height of given tree 
reflectance of a given tree 
expected value of tree reflectance 
expected value of tree height 
diameter of a given tree 
variance of reflectance 
normalized vegetation index 
mean of reflectance 
number of reflectance bands 
highest level of aggregation 
control parameter for shadow computation 
(0 = no shadows computed, 1 = compute shadows) 
control parameter for elevation computation 
(0 = compute elevations, 1 = read elevations) 
control parameter for vegetation reflection 
computation (0 = compute, 1 = read from file) 
control parameter for special case of tree 


resolution 

finest 
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VREF 


C 

C 

C 


diameter = 1 pixel (0 = dia not equal to 
one pixel, 1 = dia equal to 1 pixel) 
vegetation reflectance distribution 


C Last revised: 1/7/88 


REAL XI (20000) ,X2(20000) ,RNODE(2, 150, 150) ,TDIA(20000) 

TREF(2, 20000) , ELEV (-10 : 150, 150), ' 

LAM, ETDIA, ETREF (2) , STREF (2) , THGT (20000) , 

MREF (2) , RAGNODE (2, 150, 150) , 

DIST_SQ,DIST (150) , 

PCTV (150, 150) , PCTVAG (150, 150) , ALPHA, THETA, PCTS (150 1501 

PCTSAG (150, 150) , PCTVI (150, 150) , PCTviAG (150, 150) ' 

PCTGS (150,150) , PCTGSAG (150, 150) ,NVI (150, 150) SI* S2 
PCTVS (150,150) ,PCTVSAG (150, 150) , VREF (150, 150) , VI, V2 

INTEGER NUM, IDIM, IBAND, IA, IAG, IELEV, ISHAD, TVREF, IVDIA 

OPEN (10 / FILE= , CAN14 . IN ? , STATUS= 1 0LD 1 ) 

FILE= ’ CANl 4 . OUT • , STATUS= ' NEW 1 , FORM= * FORMATTED • ) 

OPEN (13, FILE- ’ Cl 4BPL . DAT ’ , STATUS= ’ NEW 1 , CARRIAGECONTROL= ’LIST') 

C Input model parameters 


& 

& 

& 

& 

& 

& 

& 

& 


READ (10,*) 
READ (10,*) 
READ(10,*) 
READ (10, *) 
READ(10,*) 
READ (10,*) 
READ(10,*) 
READ(10,*) 
READ (10,*) 
READ (10,*) 
READ (10,*). 


IBAND 
IDIM, IAG 
IELEV, ISHAD 
IVREF, IVDIA 
VI, V2 
SI, S2 
LAM 

ALPHA, THETA 
ETHGT, STHGT 
ETREF(l) , STREF (1) 
ETREF (2) , STREF (2) 


C Call RANGEN to generate tree locations, heights, and diameters 

CALL RANGEN (IELEV, IDIM, LAM, ALPHA, ETHGT, STHGT XI X2 
& THGT,TDIA, NUM, IBAND, ELEV) ' ' ' 


C Call OUTPAR to write inputted parameters and number of trees. 

CALL OUTPAR (LAM, ALPHA, THETA, ETHGT, STHGT, ETREF, STREF, IBAND, NUM) 
C Call ELEVA to compute elevation of tree at each node. 


CALL ELEVA (ISHAD, IELEV, IDIM, NUM, ELEV, ALPHA, XI, X2, THGT, TD IA) 
C Call OUTEL to print elevation at each node. 


CALL OUTEL (IDIM, ELEV) 


r cfi REFL f C t0 com P ute the reflectance at each node (before 
C shadow subroutine) . ' s 

CALL REFLEC (IDIM, NUM, XI, X2, TDIA, ETREF, STREF, SI , S2 , 

& TREF , RNODE , IBAND , PCTV, IVREF , IVDIA, VREF , VI , V2 ) 


223 


noon 


C Call SHADOW to compute shadowed pixels (reflectivity = 0.0). 

CALL SHADOW (ISHAD, IDIM, ALPHA, THETA, ETHGT, STHGT, 

& RNODE, ELEV, IBAND, PCTS, PCTV, PCTVI, PCTGS) 

C Compute average reflectance for increasing levels of aggregation 

DO 100 IAGR =1,4 
IA = 1 

IF (IAGR.EQ.2) THEN 
IA = 5 

ENDIF 

IF (IAGR. EQ. 3) THEN 
IA = 10 
ENDIF 

IF (IAGR.EQ. 4) THEN 
IA = 30 
ENDIF 

IF (IA.LT.2) THEN 

CALL OUTREF (IDIM, NUM, RNODE, IA, 

5 IBAND, PCTV, PCTS) 

CALL OUTPLOT (IDIM, RNODE, IA, IAG, IBAND, PCTV, PCTVI, 

6 PCTVS, PCTS, PCTGS, ELEV) 

ELSE 

CALL AGGREG (IDIM, RNODE, RAGNODE, IA, IBAND, PCTV, PCTVAG, 

& PCTS, PCTSAG, PCTVI, PCTVIAG, PCTGS, PCTGSAG) 

C call OUTREF to print array of reflectances and subpixel 
C components at each node. 

CALL OUTREF (IDIM, NUM, RAGNODE, IA, 

& IBAND , PCTVAG, PCTSAG) 

CALL OUTPLOT ( IDIM, RAGNODE , IA, IAG, IBAND , PCTVAG, 

& PCTVIAG, PCTVSAG, PCTSAG, PCTGSAG, ELEV) 

ENDIF 

C Subroutine to summarize statistics of subpixel percentages 

CALL OUTPCT (IDIM, IA, RAGNODE, PCTVAG, PCTVIAG, PCTSAG, PCTGSAG) 

100 CONTINUE 
STOP 
END 


This subroutine aggregates the reflectance values and subpixel 
percentages up to the level IA 

SUBROUTINE AGGREG (IDIM, RNODE, RAGNODE, IA, IBAND, PCTV, PCTVAG, 

& PCTS, PCTSAG, PCTVI, PCTVIAG, PCTGS, PCTGSAG) 

REAL RNODE (2, 150,150) , RAGNODE (2, 150, 150) , PCTV (150, 150) , 

& PCTVAG (150, 150) , PCTS (150, 150) , PCTSAG (150, 150) , PCTVI (150, 150) , 
& PCTVIAG (150, 150 ), PCTGS (150, 150) , PCTGSAG (150, 150) 

INTEGER IDIM, IA, IBAND, IX, IY 

C Initialize variables 
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IX * 0 
IY - 0 

DO 50 I = 1, IDIM 
DO 50 J » 1 , IDIM 

DO 60 IB = 1, IBAND 

RAGNODE (IB, I, J) = 0. 
60 CONTINUE 

PCTVAG (I, J) = 0. 

PCTSAG (I, J) = 0. 

PCTVIAG (I, J) =0.0 
PCTGSAG (I, J) =0.0 
50 CONTINUE 


C Compute aggregated matrix 


DO 100 I = 1, IDIM-IA+1, IA 
IX = IX + 1 

DO 110 J = 1, IDIM-IA+1, IA 
IY = IY + 1 
DO 120 II = 1, IA 
DO 120 J1 = 1,LA 

DO 130 IB = 1, IBAND 

RAGNODE (IB, IX, IY) = RAGNODE (IB, IX, IY) + 


130 CONTINUE 

PCTVAG (IX, IY) = 

& 

PCTSAG (IX, IY) = 

& 

PCTVIAG (IX, IY) = 

& 

PCTGSAG (IX, IY) = 


120 CONTINUE 

110 CONTINUE 
IY = 0 

100 CONTINUE 
RETURN 
END 


RNODE (IB, I+Il-l, J+Jl-1) / (LA**2) 

PCTVAG (IX, I Y) + 

PCTV (I+Il-l, J+Jl-1 ) / (IA**2) 

PCTSAG (IX, IY) + 

PCTS (I+Il-l, J+Jl-1) / (IA**2) 

PCTVIAG (IX, IY) + 

PCTVI (I+Il-l, J+Jl-1) / (IA**2) 
PCTGSAG (IX, IY) + 

PCTGS(I+I1-1, J+Jl-1)/ (IA**2) 


Subroutine RANGEN samples from appropriate random 
distributions to generate tree center locations, 
then tree heights and diameters. 


Variables: 

RAD 

NUM 

SECNDS 

GGUBFS 

GGNML 


radius of tree 
number of trees 

MICROVAX function returns time in sec. 
IMSL uniform random number generator 
IMS1 normal random number generator 


& 


SUBROUTINE RANGEN (IELEV, IDIM, LAM, ALPHA, ETHGT, STHGT 

XI , X2 , THGT, TDIA, NUM, IBAND, ELEV) 
REAL LAM,ETHGT, STHGT, MDIA, ELEV (-10 : 150, 150) , 
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& RADIUS, OMEGA, TDIA(20000),THGT(20000) , 

& XI (20000), X2 (20000) ,MHGT, 

& GGUBFS,SECNDS,R(1) , OVERLAP 

INTEGER NR, IBAND, 1X1 (20000) , 1X2 (20000) 

DOUBLE PRECISION SEED 
SAVE SEED 

PARAMETER (PI=3 . 14159) 

C Zero variables. 


RADIUS =0.0 
NUM = 0 
NR = 1 


C 

C 


If(IELEV=l) then read elevations computed outside 
this program. 

IF ( IELEV . GT . 0 . 5 ) THEN 
READ (1, *) 

DO 200 J=1 , IDIM 


200 


READ (1, *) 
READ(1, *) 
READ (1, *) 
READ (1, *) 
READ(1,*) 
READ (1, *) 
READ (1, *) 
CONTINUE 


(ELEV (I, J) , 1= -10,0) 
(ELEV (I, J) , 1= 1,25) 

(ELEV (I, J) , 1= 26,50) 

(ELEV (I, J) , 1= 51,75) 

(ELEV (I, J) , 1= 76,100) 

(ELEV (I, J) , 1= 101,125) 
(ELEV (I, J) , 1= 126,150) 


DO 210 J=l, IDIM 

DO 210 I=-10, IDIM 

IF (ELEV (I, J) .GT. 0 . 0) THEN 
NUM = NUM + 1 
THGT(NUM) = ELEV (I, J) 
XI (NUM) = I 
X2 (NUM) = J 


TDIA(NUM) =1.0 
ENDIF 

210 CONTINUE 
GO TO 110 
ENDIF 

C Obtain seed from internal clock. 

SEED = SECNDS (0.0) + 10000. 

C Generate tree centers 

C Simulate 2-D Poisson Process 
C Source: Cox, Point Processes 

C Method: generate variable X = PI* (R2**2 - Rl**2) from 

exponential distribution, and angle, 

OMEGA, from uniform distribution on 0,2PI 

To obtain exponential deviate, X, from U uniform (0,1) * 
X = (-1/R) *LOG (U) ; 
where f(x) = R * exp(R*x) 

GGUBFS(SEED) - IMSL uniform random number generator 
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c 


Maximum radius of tree location equal to extent 
plus ten times tree diameter. 


of field 


C Maximum radius is SQRT (2) *IDIM/2 + 3*RMAX 


C 


RADMAX 
10 RADIUS 
& 


3 . * (ETHGT* (TAN (ALPHA*6 . 283/360) ) *2 . ) + 106 1 
SQRT ( (- (1 . /LAM) *LOG (GGUBFS (SEED) ) ) /PI 

+ RADIUS*RADIUS) 


If radius exceeds extent of field return to main program 


IF (RADIUS. GT. RADMAX) GO TO 20 
IF (NUM.GE. 20000) GO TO 20 


C Increment tree number, generate angular coordinate. 


NUM = NUM+1 

OMEGA = GGUBFS (SEED) * 2. * PI 


1X1 (NUM) = RADIUS * COS (OMEGA) + 75 
1X2 (NUM) = RADIUS * SIN (OMEGA) + 75 ] 

C Put tree at center of pixel by truncating to integer value. 

XI (NUM) = IXl(NUM) 

X2 (NUM) = 1X2 (NUM) 


C 

C 


Generate tree height (normal distribution) . . 

Set bounds for acceptable range of heights (l-io ‘meters ’only) .* 

CALL GGNML(SEED,NR,R) 

THGT (NUM) = ETHGT + (STHGT*R(1) ) 

IF (THGT (NUM) .LT. 0 . ) THEN 
THGT (NUM) =1.0 
END IF 

IF (THGT (NUM) .GT. 10.0) THEN 
THGT (NUM) =9.9 
ENDIF 

MHGT = MHGT + THGT (NUM) 


C Compute tree diameter 

TDIA(NUM) = THGT (NUM) *TAN( ALPHA* 6. 283/360) *2 
MDIA = MDIA + TDIA(NUM) 


C 

c 


Compute inhibitory field, 
then no tree is generated 


II tree locations overlap 
(Decrement tree NUM) . 


& 


IMAX = MIN (NUM, 500) 
ITEST = 0. 

IF (IMAX. LT. 2) GO TO 600 
DO 600 IT = 1, IMAX-1 


xr UXE5T.GT.0.1) GO TO 600 
DIST_TR = SQRT ( (XI (NUM) - XI (NUM- IT) ) * 
(X2 (NUM) - X2 (NUM-IT) ) ** 2 ) 
OVERLAP = ( (TDIA (NUM) +TDIA (NUM-IT) )/2) 
IF (OVERLAP. GT. 0.1) THEN 
NUM = NUM-1 


+ 

DIST TR 
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ITEST -1.0 
ENDIF 

600 CONTINUE 

IF (ITEST. GT. 0. 1) GO TO 10 

GO TO 10 

20 MHGT = MHGT/NUM 
MDIA = MDIA/NUM 
110 RETURN 
END 


Subroutine ELEVA calculates the elevation at each pixel based 
on the tree center height and cone angle. Ground surface is 
assumed flat and horizontal at elev = 0.0. 

SUBROUTINE ELEVA ( ISHAD, IELEV, IDIM, NUM, ELEV, 

& ALPHA, XI, X2,THGT,TDIA) 

REAL ELTEMP,ELEV (-10:150, 150) , ALPHA, XI (20000) ,X2 (20000) , 

& THGT (20000) ,TDIA (20000) , XX, YY 

INTEGER IDIM, NUM, ITDIA,X11,X22 

C If (IELEV-1) then use elevations computed outside this program 
C and skip this subroutine 

IF (IELEV. GT. 0.5) GO TO 110 

C Initialize pixels at 0.0 elevation. 

DO 20 I- -10, IDIM 
DO 20 J- 1, IDIM 
ELEV (I, J) =0.0 
20 CONTINUE 

C If (ISHAD-0) then above initializations stand and skip 
C the rest of this subroutine 

IF (ISHAD.LT. 0.5) GO TO 110 

TANCON = TAN (ALPHA*6. 283/360) 

C Iteration to compute elevation 

DO 100 N=1,NUM 
ITDIA = TDIA(N) 

DO 50 I—ITDIA, ITDIA 
DO 50 J—ITDIA, ITDIA 

C Find nearest integer pixel by truncating 

Xll = XI (N) + I 
X22 = X2(N) + J 

IF (Xll. GT. IDIM .OR. X11.LT.-10) GO TO 50 
IF (X22.GT.IDIM .OR. X22.LT.1) GO TO 50 
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C Compute elevation based on distance from tree center 

ELTEMP = THGT (N) - SQRT ( ( (Xll-Xl (N) ) **2) + 
& ( (X22-X2 (N) ) **2) ) /TANCON 

C For overlapping tree canopies, take highest canopy 

IF (ELTEMP. LT.ELEV (XI 1,X22) ) THEN 
GO TO 50 
ELSE 

ELEV (Xll, X22 ) = ELTEMP 
ENDIF 

50 CONTINUE 

100 CONTINUE 
110 RETURN 
END 


C 

C 

C 


Subroutine OUTEL prints the elevation of each pixel 
SUBROUTINE OUTEL (ID IM, ELEV) 


800 

850 

870 


REAL ELEV (-10: 150, 150) 

INTEGER IDIM 
WRITE (11,800) 

FORMAT (/ / , 2X, 'Pixel elevation in meters',/) 
DO 850 J=l, 30 

WRITE (11, 870) (ELEV(I,J) ,1=1,15) 

CONTINUE 
FORMAT (15F6. 2) 

RETURN 

END 


Subroutine REFLEC computes reflectance at each node based 
on superposition of trees on soil background. 

Variables : 

DIST_SQ - square of distance from node to tree center 


SUBROUTINE REFLEC (IDIM, NUM, XI, X2, TDIA, ETREF, STREF, SI S2 
& TREF, RNODE, IBAND, PCTV, TVREF, IVDIA, VREF^ VI ' V2 ) 

REAL XI (20000) , X2 (20000) , TDIA (20000) , ETREF (2) , STREF (2 ) 

& TREF (2, 20000) , RNODE (2, 150, 150) ,DIST SQ,PCTV(150, 150) 

4 SECNDS, R (1 ) ,MREF (2 ) , VREF (150,150) , vT, V2, SI, S2 

INTEGER IDIM, NUM, INUM, NR, IVREF, IVDIA, 1X1, 1X2 
DOUBLE PRECISION SEED 
SAVE SEED 


SEED = SECNDS (0.0) + 10000. 
NR = 1 


C Initialize percent vegetation 

DO. 15 1 = 1, IDIM 
DO 15 J = 1, IDIM 
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PCTV (I, J) = 0. 

15 CONTINUE 

C Soil distribution obained from the TURNING BANDS model 
C computed outside this program. Just read input here. 

DO 40 1=1, IDIM 

READ (2, *) (RNODE ( 1 , I , J) , J= 1,15) 

READ (2,*) (RNODE (1,1, J) , J= 16,30) 

READ (2,*) (RNODE (1,1, J),J= 31,45) 

READ (2,*) (RNODE (1, 1, J) , J= 46,60) 

READ (2,*) (RNODE (1,1, J) , J= 61,75) 

READ (2,*) (RNODE (1,1, J) , J= 76,90) 

READ (2,*) (RNODE (1,1, J) , J= 91,105) 

READ (2, * ) (RNODE (1, I, J) , J=106, 120) 

READ (2,*) (RNODE (1, 1, J) , J=121, 135) 

READ (2,*) (RNODE (1, 1, J) , J=136, 150) 

C Infrared soil reflectance linearly related to visible reflectance 
DO 40 J=l, IDIM 

IF (RNODE (1, I, J) .LT. 0.05) THEN 
RNODE (1, I, J) = 0.05 
ENDIF 

RNODE (2, I, J) = (Sl*RNODE (1, I, J) ) + S2 
40 CONTINUE 

C If ( IVREF . EQ . 1 ) then read visible tree reflectance distribution 
C from outside this program and assign reflectance at tree centers 

IF ( IVREF . GT .0.5) THEN 
DO 200 1=1, IDIM 

READ (3,*) (VREF (I, J) , J= 1,15) 

READ (3,*) (VREF (I, J) , J= 16,30) 

READ (3,*) (VREF (I, J) , J- 31,45) 

READ (3, *) (VREF (I, J) , J= 46,60) 

READ (3, *) (VREF (I, J) , J= 61,75) 

READ (3, *) (VREF (I, J) , J= 76,90) 

READ (3, *) (VREF (I, J) , J= 91,105) 

READ (3,*) (VREF (I, J) , J=106, 120) 

READ (3,*) (VREF (I, J) , J=121, 135) 

READ (3, *) (VREF (I, J) , J=136, 150) 

200 CONTINUE 

DO 250 IN=1 , NUM 
1X1 = XI (IN) 

C (can't handle negatives in some arrays) 

IF (1X1 . LT .1.0) GO TO 250 
1X2 = X2 (IN) 

TREF (1, IN) = VREF (1X1, 1X2) 

C Assume IR tree reflectance logrithmicly related to VIS reflectance 

TREF (2, IN) = (V1*TREF (1, IN) ) + V2 

250 CONTINUE 
GO TO 270 
ENDIF 
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c 

c 

c 


60 


«cfS«!'Rlvu[cn%^/2i r / 8 f 7 e ^ anCe ‘? 0mal attribution, for 
same random relation) between visibTe^and^nf rel J tlon (therefore 
Therefore only one random 

DO 60 IB = 1, IBAND 
DO 60 INUM = 1, nuM 
CALL GGNML (SEED, NR, R) 

IB = 1 

TREF (IB, INUM) = ETREF (IB) + (STREF (IB) *R ( 1 ) ) 

®EF(IB) a MBEF (IB, * TF£F(IB,lS ’ R(1)) 

IB — 2 

TREF (IB, INUM) = ETREF (IB) + (STREF (IB) *R < 1 n 
C0 S® - * TBEFdB.lS 10 ' R<1,) 


270 CONTINUE 


C 

C 

C 




C equal^tre^reflectanS^nlv a^tha^' ^ 6n pixel reflec tance 

C time-consuming calculation 17 * 9rid; Skip followin g 


290 

295 


IF (IVDIA. GT. 0 . 5) THEN 
DO 295 1=1, NUM 


1X1 = XI (I) 

IF (1X1 .LT.l . 0) GO TO 295 
1X2 = X2 (I) 


DO 290 IB=1, IBAND 
RNODE(IB,IXl,IX2)= 
CONTINUE 


TREF (IB, I) 


PCTV (1X1, 1X2) =1.0 
CONTINUE 
GO TO 280 
ENDIF 


DO 10 I = 1, NUM 
Xll = XI (I) 

X22 = X2 (I) 

DO 20 J = 1, IDIM 

DO 20 K = 1, IDIM 
Z1 = J 
Z2 = K 


C 

C 


DIST_SQ - (Xll - Zl)**2 + (X22 - 
IF ( (TDIA ( I ) / 2 ) . GE . SQRT (DIST_SQ) ) 


Z2) **2 
THEN 


of'h S W n S t SUnied Phat reflectan ce of the overlap portion 
of two trees is equal to the latter generate/tree. 


50 


DO 50 IB=1, IBAND 
RNODE (IB, J, K) = 
CONTINUE 
PCTV ( J, K) =1.0 
ENDIF 


TREF (IB, I) 
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20 CONTINUE 

10 CONTINUE 
280 RETURN 
END 


Subroutine to compute shadows cast by trees. Sun is 
assumed directly from west (left of page) . 

SUBROUTINE SHADOW (ISHAD, IDIM, ALPHA, THETA, ETHGT, STHGT, 

& RNODE , ELEV, IBAND , PCTS , PCTV, PCTVI , PCTGS ) 

REAL ALPHA, THETA, ETHGT, STHGT, RNODE (2, 150, 150) , 

& ELEV (-10 : 150, 150) , TANCON, TANSUN, PCTS (150, 150) , 

& PCTV (150, 150) , PCTVI (150, 150) , PCTGS (150, 150) 

INTEGER ISHAD, IMAX, IBAND 

C Initialize variables. 

DO 100 1=1, IDIM 
DO 100 J=l, IDIM 
PCTS (I, J) = 0. 

PCTVI (I, J) = PCTV (I, J) 

PCTGS (I, J) = 0. 

100 CONTINUE 

C If there are no shadows (ISHAD = 0.0), above initializations 
C stand and rest of algorithm can be skipped. 

IF (ISHAD.LT. 0.5) GO TO 60 

TANCON = TAN (ALPHA*6. 283/360) 

TANSUN = TAN ( (90-THETA) *6.283/360) 

IMAX = (ETHGT + (5*STHGT) ) /TANSUN 

DO 50 J = 1 , IDIM 

DO 50 I = IDIM, 1,-1 
DO 40 IN = 1, IMAX 

IF ( (I-IN) .LT.-10) GO TO 40 

IF ( ( ( (ELEV ( (I-IN) , J) ) - ELEV (I, J) ) /IN) .GE. TANSUN) THEN 
RNODE (1, I, J) =0.0 
RNODE (2, I, J) =0.0 

C Initialize percent shadow 

PCTS (I, J) =1.0 
IF (PCTV (I, J) .GE. 0.95) THEN 
PCTVI (I, J) =0.0 
ELSE 

PCTGS (I, J) =1.0 
ENDIF 



ENDIF 

40 

CONTINUE 

50 

CONTINUE 

60 

RETURN 


END 
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This subroutine prints the reflectances 
and a summary of the statistics 


at each node 


SUBROUTINE OUTREF (IDIM, NUM, RNODE, IA, 

IBAND , PCTV, PCTS ) 

REAL RNODE (2, 150, 150) , 

LAG, DIST, PCTV (150,150) , PCTS (150, 150) 
INTEGER IDIM, NUM, IBAND, IA, ITAU, IMAX, IMIN 
ITAU = IDIM/IA 
IMAX - 15 


IMIN = MIN (IMAX, ITAU) 

WRITE (11, 805) IA 

805 FORMAT (//,' 

& 

& '/, 26X, 'Level of aggregation =' 14) 

DO 810 IB = 1, IBAND 
WRITE (11, 815) IB 

815 STJ j<?MN lxel leflectance: 

830 Co™' 11,850 ’ (RN0DE(IB ' I - J >- I-l.IMIN) 

850 FORMAT (15F6. 2) 

810 CONTINUE 

WRITE (11, 845) ’Percent vegetation' 

845 FORMAT (//,2X, A,/) 

DO 840 J =1, IMIN 

Qan ^ X 3 ITE(11 ' 850 ) (PCTV ( I, J), 1=1, imin) 

840 CONTINUE 

WRITE (11, 845) 'Percent shadow' 

DO 842 J =1, IMIN 

"*^<11,850) (PCTS ( I, J), 1=1, IMIN) 

842 CONTINUE 


RETURN 

END 


C 

c 

c 

c 


This subroutine writes inputted storm 
number of trees generated. 


parameters and 


SUBROUTINE OUTPAR (LAM, ALPHA, THETA, ETHGT, STHGT, 

& „„ ETREF,STREF, IBAND, NUM) 

REAL LAM, ALPHA, THETA, ETHGT, STHGT, ETREF (2) , STREF (2) 
INTEGER IBAND, NUM 11 

WRITE (11, 860) 'CANOPY MODEL: VERSION 11', 

* 'LAMB', 'ALPH', 'THET', ' E [THGT] ' , ' S [THGT1 ' 

0 ,. S , ’E(REF_VIS] ’, 'S[REF VIS] ', 'EfREF IR1 ' ' S fRE 

860 FORMAT (//, 2X, A, / , 2X, A, 2X, A, 2X, A~2X, A, 2X, A, 1X~A, 2X*A 2X^ 
WRITE (11, 870) LAM, ALPHA, THETA, ETHGT, STHGT X ' A '^ X ' A ' 2X ' 
& (ETREF (IB) , STREF (IB) , IB-1, IBAND) 

870 FORMAT (/,F6.3,3F6.1,F10.3,F11.3,F12.3,F10.3,F11 3) 

WRITE (11, 820 ) NUM ' 

820 FORMAT (//, 2X, ’Number of trees =',X 18) 

RETURN ' 

END 


IR] ’ 
2X, A) 
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This subroutine computes subpixel components 
and arranges output for use in the DRAW package 
of the microVAX II. 

SUBROUTINE OUTPLOT (IDIM, RNODE, IA, I AG, 

& IBAND , PCTV, PCTVI , PCTVS , PCTS , PCTGS , ELEV ) 

REAL RNODE (2, 150, 150) , LAG, DIST, PCTV (150, 150) , 

& PCTS (150, 150), PCTVI (150, 150), PCTGS (150, 150) , 

& PCTG (150, 150) , PCTGI (150, 150) , PCTVS (150, 150) , 

& ELEV (-10: 150, 150), NVI (150, 150) 

INTEGER IDIM, IA, LAG, ITAU, IMAX, IMIN, IPLOTB 
ITAU = I DIM/ 1 A 
IMAX = 15 

IMIN = MIN (IMAX, ITAU) 

C Initialize variables 

DO 50 J = 1 , ITAU 

DO 50 I = 1, ITAU 
PCTVS (I, J) => 0. 

PCTG (I, J) = 0. 

PCTGI (I, J) = 0. 

NVT (I, J) = 0. 

50 CONTINUE 

C Write header for draw files (once) 

IF (IA.GT.l) GO TO 960 
IPLOTB =12 
WRITE (13, 970) IPLOTB 
970 FORMAT (' CANOPY MODEL: VERSION 14', 

& /, 14, /, ' VIS_REF ' , / , ' INF_REF ' , / , 'M',/, 'MI 1 , /, 'MS * , /, 'MI/M', 
& /, 'MS/M',/, 'G',/, 'GI',/, 'GS',/, 'S',/, ’NVI') 

C Skip printing for aggregations less than 3 

960 IF (IA.LT.3) GO TO 980 

C Print file or reflectance and percent vegetation 
C for last aggregation. 

WRITE (13, 700) IA 

700 FORMAT (/, ' Level of aggregation is', 14,/) 

DO 100 J=1 , ITAU 
DO 100 1=1, ITAU 

PCTVS (I, J) = PCTS (I, J) - PCTGS (I, J) 

PCTVI (I, J) = PCTV (I, J) - PCTVS (I, J) 

PCTG (I, J) = 1.0 - PCTV (I, J) 

PCTGI (I, J) = PCTG (I, J) - PCTGS (I, J) 

RSUM = RNODE (2, I, J) + RNODE (1, I, J) 

IF (RSUM.EQ. 0 . 0) THEN 
NVI ( I , J) = 0.0 
GO TO 150 
ENDIF 
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NVI(I,J) = (RNODE (2, I, J) - RNODE (1, 1, J) ) /RSUM 

C Correction for zero-divide (at low aggregations) 

150 IF (PCTV(I,J) .LT. 0.0001) THEN 

PCTV (I, J) = .0001 
ENDIF 


WRITE (13, 200) 

& 

& 

& 

& 

& 

100 CONTINUE 
200 FORMAT (12F8. 3) 


100*RNODE (1,1, J) , 100*RNODE (2, 1, J) 
100*PCTV (I, J) , 100*PCTVI (I, J) , 

100*PCTVS (I, J) , 100*PCTVI (I, J) /PCTV (I j) 
100*PCTVS (I, J) /PCTV (I, J) , 100*PCTG (I J) 
100*PCTGI (I, J) , 100*PCTGS (I, J) , 

100*PCTS (I, J) , 100*NVI (I, J) 


980 RETURN 
END 


C 

c 

c 


Subroutine t0 suinmari2 e statistics of the subpixel percentages' 
SUBROUTINE OUTPCT (IDIM, IA, RNODE, PCTV, PCTVT, PCTS, PCTGS) 

REAL PV, PVI, PVS, PG, PGS, PGI, PS, PCTV (150, 150) , PCTS (150 150) 

‘ 50 - 150 > ■ — ^ ■ 


IMAX = 15 

ITAU = IDIM/IA 

IMIN = MIN (IMAX, ITAU) 


PV - 0. 
PVI = 0. 
PGS = 0. 
R1 - 0. 
R2 = 0. 
PS = 0. 


100 


DO 100 1=1, ITAU 
DO 100 J=l, ITAU 

PV = PV + PCTV (I, J) / (ITAU**2) 

PVI = PVI + PCTVT (I, J) / (ITAU**2) 

PGS = PGS + PCTGS (I, J) / (ITAU**2) 

R1 = R1 + RNODE (1,1, J) / (ITAU**2) 

R2 = R2 + RNODE (2, I, J) / (ITAU**2) 

PS = PS + PCTS (I, J)/(ITAU**2) 

IF ( (RNODE (2,1, J) + RNODE (1,1, J) ) .EQ 0 0) 

NNVI =0.0 
GO TO 100 
ENDIF 

NNVI = NNVI + ( (RNODE (2, 1, J) -RNODE (1, 1, J) ) / 

CONTINUE (RNODE (2, 1, J) +RNODE (1, 1, J) ) )/(ITAU**2) 


THEN 


C Correction for zero divide 
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IF (PV.LT. 0.0001) THEN 
PV =» 0.0001 
ENDIF 

WRITE (11, 210) ' Subpixel component averages 'R1 'R2 ' , 

& 'M' , 'MI' , 'MS' , ' MI/M' , 'MS/M', ' G ' , ' GI ' , ' GS ' , 'S', 1 NVI ' 

210 FORMAT (//,' 

& /, A , /, 4X, A, 4X, A, 5X, A, 4X, A, 4X, 

&A, 2X, A, 2X, A, 5X, A, 4X, A, 4X, A, 5X, A, 4X, A, //) 


PG = 1-PV 
PGI = PG - PGS 

WRITE (11, 200) R1 , R2 , PV, PVI , PV-PVI , PVT /PV, (PV-PVI) /PV, 
& PG, PGI, PGS, PS,NNVI 

200 FORMAT (2F6.3, 9F6.3,F7.2) 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


INVERSE PROCEDURE FOR CASE II 
(Program M_II.FOR) 

SfSo”a C ^?o u tp“1ne d ?or'a m ^cS?eS^'° f % 9r ° Up ° f 
in tha red-lR scattergrain, identified^by^slop^A^interceo? 9 f 

and Bandwidth BANDWIDTH. Program computes ? ' intercept ' B ' 

a ^ lmate of P erce nt vegetation cover, pm EST from simulan 
and compares it to mean of actual percent " S 


REAL RHODE < 2 ) , PM ACT, BIGIR, SMALLIR, BIGVS, SMALLVS 
SMALLPM, BIGPM, ' 

RIR_SUM, MEANIR, SDEV VSSL, SDEV IRSL 
RIRSQ_SUM, RVS_SUM, MEANVS , RVSSQ SUM," 

BANDWIDTH, A, B,CVP ARM (10 ) , ~ 

PM_EST_VS , PM_EST_IR, PARM, MACT ( 300 ) 

INTEGER NN,NUM^^^ ,PM - SDEV ' PM - SUM ' PM - SSQ ' NUM ' 
OPEN (10, f ILE='m_II. IN’ , STATUS= ' OLD ' ) 

ASSIGN INPUT DRAW FILE TO FOROll 
OPEN(12,FILE='M_II .OUT' , STATUS='NEW’ ) 

OPEN ( 1 3 , FILE= ' M_I IL . OUT ' , STATUS* ' NEW • ) 


SMALLIR = 500.0 
SMALLVS = 500.0 
SMALLPM = 500.0 
NUM =0.0 
RIR_SUM =0.0 
RIRSQ_SUM =0.0 
RVS_SUM =0.0 
RVSSQ_SUM =0.0 
PM_SSQ = 0.0 
PM_SUM =0.0 
NN = 0 


C 

c 

c 

c 


c 

c 


Raad input parameters 

RA a H R fJ D(10 ^* ) N^AIRS^B, BANDWIDTH, PARM 
Read sdev of soil line 

READ (10,*) SDEV_VSSL, SDEV IRSL 
Skip header ~ 

READ (11, 300) IPLOT 
300 FORMAT (7,14,////////////) 

Iterate for each data pair 
DO 100 I = 1, NUMPAIRS 


READ (11,*) (RNODE (K) , K= 
PM_ACT = CVP ARM (PARM) 
Compute perpendicular distance 
the line RIR = A*RVIS + B 


1,2), ( CVP ARM (J) , J=l, 10) 

between a chosen point and 


C 


AP - TAN (- (3. 142/2 - ATAN (A) ) ) 
BP - RNODE ( 2 ) - (AP*RNODE (1) ) 
RVS = (B - BP) / (AP - A) 

RIR = A* RVS + B 
DIS =SQRT ( (RIR-RNODE (2) ) **2 + 
WRITE (5,*) AP, BP, RIR, RVS, DIS 


(RVS 


RNODE (1) ) **2) 


IF (DIS . GT. BANDWIDTH) GO TO 100 
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NUM = NUM + 1.0 


NN - NN +1 
MACT(NN) =» PM_ACT 

C Sum IR data 

RIR_SUM = RIR_SUM + RNODE (2) 
RIRSQ_SUM = RIRSQ_SUM + (RNODE (2) **2) 
BIGIR = MAX (BIGIR, RNODE (2) ) 

SMALLIR = MIN ( SMALLIR, RNODE ( 2 ) ) 

C Sum Visible data 


RVS_SUM = RVS_SUM + RNODE (1) 

RVSSQ_SUM = RVSSQ_SUM + (RNODE (1) **2) 

BIGVS * MAX (BIGVS, RNODE (1) ) 

SMALLVS =* MIN ( SMALLVS , RNODE ( 1 ) ) 

C Sum percent canopy cover 

PM_SUM = PM_SUM + PM_ACT 

PM_SSQ = PM_SSQ + PM_ACT**2 

BIGPM = MAX (BIGPM, PM_ACT) 

SMALLPM = MIN (SMALLPM, PM_ACT) 

100 CONTINUE 

MEANIR = RIR_SUM/NUM 

SDEVIR = SQRT ( (RIRSQ_SUM - NUM* (MEANIR**2) ) / (NUM-1) ) 
MEANVS = RVS_SUM/NUM 

SDEWS = SQRT ( (RVSSQ_SUM - NUM* (MEANVS**2) ) / (NUM-1) ) 
PM_MEAN = PM_SUM/NUM 

PM_SDEV = SQRT ( (PM_SSQ - NUM* (PM_MEAN**2) ) / (NUM-1) ) 

C Write input and output parameters 

WRITE (12, 210) 'NUM', 'A', 'B', 'BANDWIDTH' , NUM, A, B, BANDWIDTH 
WRITE (5, 210) 'NUM', 'A', 'B', ' BANDWIDTH' , NUM, A, B, BANDWIDTH 
210 FORMAT (5X, A, 7X, A, IX, A, 2X, A, / ,F8 .0, 2F8 . 3, F10 . 3) 

C Write summary of reflectance statistics 

WRITE (12,200) 'BAND' , 'MEAN' , 'SDEV', 'SD/MN', 'MAX' , 'MIN' , 

& ' IR' , MEANIR, SDEVIR, SDEVIR/MEANIR, BIGIR, SMALLIR, 

& 'VIS' , MEANVS, SDEWS, SDEWS/MEANVS, BIGVS, SMALLVS 

WRITE (5, 200) 'BAND', 'MEAN', 'SDEV', 'SD/MN', 'MAX', 'MIN', 

& ' IR' , MEANIR, SDEVIR, SDEVIR/MEANIR, BIGIR, SMALLIR, 

& 'VIS' , MEANVS, SDEWS, SDEWS/MEANVS, BIGVS, SMALLVS 

200 FORMAT (/, 4X,A, 4X, A, 4X, A, 3X, A, 5X, A, 5X, A, //, 

& 6X,A, 5F8.2,/,5X,A, 5F8.2) 

C Compute canopy cover using variances 
C COMPUTE E{GI} 


PM_EST_VS = 100* (1 - ( SDEWS / SDEV_VSSL ) ) 
PM EST IR = 100* (1 - ( SDEVIR/ SDEV IRSL) ) 


C Write summary of canopy cover statistics 

WRITE (*,301) ' M_ACT ' , 'SDEV', 'M_EST(VIS) ', *M_EST(IR) ', 

& PM_MEAN, PM_SDEV, PM_EST_VS, PM_EST_IR 

WRITE (12, 301) ' M_ACT ' , 'SDEV' , ’E[M]_(VIS) ', 'E[M]_(IR) ' , 

& PM_MEAN, PM_SDEV, PM_EST_VS, PM_EST_IR 

301 FORMAT ( / , 5X, A, 6X, A, 2X, A, 3X, A, / , 2F10 . 2, 2F12 . 2) 

DO 120 K=1,NN 

WRITE (13, 302) K, MACT(K), PM_EST_VS, PM_EST_IR 
120 CONTINUE 

302 FORMAT (18, 3F8. 2) 

STOP 

END 
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c 

c 

c 

c 

c 

c 


INVERSE PROCEDURE FOR CASE V, METHOD 1 
(VB2 PARV02 . FOR) 

Estimates parameters from moment analysis 
Three unknowns, EM, ETA, AT 


c 


REAL EM, EM_MIN, EM_MAX, DEM, VM, 

EGI, EGS, VGI, CMGS, CMGS_MIN, CMGS MAX, DCMGS, DCM 
AT , AT_MIN, AT_MAX , DAT ,AP, ~ 

Y1,Y2,Y3,Y(3),F,FMIN,X1,X2,X3, 

MnrJra ) m 7 ?? 1 ( 2 } ' SER ( 2 > ' SVR ( 2 > ' SCR1CR2 , RV C 2 ) , 
NDELTA,M_SAM,CMGS SAM, CMGS MIN INIT,CMGS MAX I NTT 

EM_MIN INIT, EM MAX INIT, AT~ MIN~ IN^T AT mSx^N^T 
FPRINT (10) , ATPRNT (10) , Wl, W2,W37Ai;A2 ~ ~ T ' 

ETA, ETA_MIN, ETA MAX, DETA ' ' 

INTEGER LEVEL, IDELTA 


ASSIGN FOR010.DAT MOMENT.DAT 
OPEN (11, FILE='PARVL. IN' , STATUS 
OPEN (12, FILE= 1 PARV. OUT 1 , STATUS 


— ’ OLD ' ) 
= ' NEW' ) 


C 

C 


READ (11,*) LEVEL, NDELTA, IDELTA 
READ (11,*) EM_MIN_INIT, EM MAX INIT 
READ (11,*) AT_MIN_INIT, AT MAX INIT 
READ (11,*) ETA_MIN_INIT, ETA MAX INIT 
Weighting constants for minimi zation 
READdl,*) W1,W2,W3 
Read correlation coefficients 
READ(11,*) Al, A2 
READ(11,*) A3,A4 
READ(11,*) A5,A6 
DO 800 NDAT =1,8 


C 

C 


Read fixed parameters and sample 
CALL READDAT (AP,M SAM, CMGS 
Initialize “ 


moments 
SAM, ERGI, 


VRGI , SER, SVR, SCR1R2 ) 


EMF = 100.00 
ATF = 100.00 
CMGSF = 100.00 
FMIN = 10000000000.00 
EM_MIN = EM MIN INIT 
EM_MAX = EM - MAX - INIT 
ATJMIN = AT_MIN_INIT 
AT_MAX = AT_MAX_INIT 
ETAJMIN = ETA_MIN_INIT 
ETA_MAX = ETA MAX - INIT 


460 

440 


u . U1 U 0 

CMGS_SAM = CMGS_S AM* 0. 000100 
WRITE (5, 460) ' M_S AMPLE = >,M SAM, 'C0V(M,GS) 
WRITE (12, 460) *M SAMPLE = ',M SAM, 'C0V(M GS 
FORMAT ( / / , 2 ( 3X, a7"f 9.5)) - ( ' GS 

WRITE (5, 440) 'EM' , 'AT* , 'ETA' , 'CMGS', 'Yl' 'Y 
WRITE (12, 440 ) 'EM' , 'AT' , 'ETA' 'CMGS' 'Yl' *' 
FORMAT (6X, A, 6X, A, 5X, A, 6X, A, 3 (iox , 7 X A) ' 


, CMGS_SAM 
CMGS_SAM 

’ Y3 ’ , 'FMIN' 

, *Y3', 'FMIN' 
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C Iterate over MAX and MIN values 
DO 600 IT = 1, LEVEL 

DEM = (EM_MAX - EM_MIN) /NDELTA 
DAT = (AT_MAX - AT_MIN) /NDELTA 
DETA = (ETA_MAX - ETA_MIN) /NDELTA 
EM = EM_MIN 
DO 500 IT1 = 1 , IDELTA 
EM = EM + DEM 
AT = AT_MIN 
DO 510 IT2=1, IDELTA 
AT = AT + DAT 
ATPRNT (IT2) = AT 
ETA = ETA_MIN 
DO 520 IT3=1, IDELTA 

ETA = ETA + DETA 

C Call minimization routine 

CALL FUNCT ( EM, AT , CMOS , Y , F , AP , ETA , ERGI , VRGI , 
& SER, SVR, SCR1R2 , EGI , EGS , VM, VGI , RV, W1 , W2 , W3 , 

& Al, A2, A3, A4, A5, A6, ITER) 

FPRINT (IT3) = MIN(999. 900, F) 

IF (F.LT.FMIN) THEN 
FMIN = F 
EMF = EM 
ATF = AT 
CMGSF =■ CMGS 
EGIF = EGI 
EGSF = EGS 
ETAF = ETA 
Y1 = Y (1) 

Y2 = Y (2) 

Y3 - Y (3) 

VMF = VM 
VGIF - VGI 
RV1F = RV(1) 

RV2F = RV(2) 

END IF 

520 CONTINUE 

510 CONTINUE 

500 CONTINUE 

WRITE (5,430) EMF, ATF, ETAF, CMGSF, Y1,Y2,Y3, FMIN 
WRITE (12, 430) EMF, ATF, ETAF, CMGSF, Y1,Y2,Y 3, FMIN 
430 FORMAT (3F8 . 5, F10 . 5, 3 (X, Fll . 6) ,X,F10.6) 

ITER = 0 

C Change MAX and MIN values 

EM_MAX = EMF + 0 . 4D0* (EM_MAX - EM_MIN) 

EM_MIN = EMF - 0 . 4D0* (EM_MAX - EMJMEN) 

EM_MAX = EMF + DEM 
EM_MIN = EMF - DEM 
IF (EM_MIN. LT .0.00) THEN 
EM_MIN = 0.000 
ENDIF 

IF (EM_MAX . GT. 1.00) THEN 
EM MAX = 1.000 
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ENDIF 

AT_MAX = ATF + 0 . 400* (AT__MAX - AT MIN) 
AT_MIN = ATF - 0.400* (AT_MAX - AT~MIN) 
■IF(AT_MIN.LT.0.00) THEN “ 

ATJMIN = AT MIN INIT 
ENDIF ~ “ 

ETA_MAX = ETAF + 0.400* (ETA MAX — ETA MIN) 
ETA__MIN = ETAF - 0.400* (ETA MAX - ETA~MIN) 
IF (ETA_jMIN. LT .0.00) THEN “ 

ETA_MIN = ETA MIN INIT 
ENDIF ~ ~ 

600 CONTINUE 


& 

& 

450 

Sc 


1 RV1 


”““5; 4 5°) ’EM', 'EGI 1 , 'AT', 'ETA', 'VM', 'VGI’, 'CMGS 1 
RV2 , EMF, EGIF, ATF, ETAF, VMF, VGIF, CMGSF, RV1F RV2F 
WRITE (12, 450) 'EM' , 'EGI ' , 'AT ' , ' ETA' , VM' , ' VGI" ' 04GS ' 'RV1' 
EGIF, ATF, ETAF , VMF , VGIF , CMGSF, RV1F, RV2F 
FOFMAT ,|X,A, 4X.A, |X,A. <X.A, 7X.A, SX,A, «*. A. 2 <7X,1, . /, 4 <F7 . 3, , 


800 CONTINUE 
STOP 
END 


C 


SUBROUTINE FUNCT (EM, AT, CMGS, Y, F, AP, ETA, ERGI, VRGI, 

SER, SVR, SCR1R2 , EGI , EGS , VM, VGI , RV, W1 , W2 , W3 , 

Al, A2, A3, A4, A5, A6, ITER) 

REAL Y (3) , ER (2) , VR (2) , CR1R2, EM, EGS, EGI, 

VM, VGS,VGI,W1,W2,W3,A1,A2,A3,A4,A5,A6, 

ERGI (2 ) , VR GI (2 ) , Rv {2 ) # cmgs, SER (2 ) , SVR (2 ) , SCR1R2 , 

AP , F, ETA, AT, AST, DUM, CORR, CST 


EGI = (1.00 - EM) ** (ETA + 1.00) 

EGS = 1.00 - EM - EGI 

EGSGI = 1 - EM 
X = 2.3/ ( (AP /AT) **1.50) 

VM = (X*EGSGI) + 

‘ AST - ( “^; o *^“ 0 ^< ESSGI "(-“/AP))*(1.00-X)-1.00) 

X = 2.300/ ( (AP /AST) **1.500) 

VGI = (X*EGI) + 

& (EGI** 2 . 00) * ( (EGI** (-AST/AP) ) * (1 . 00 - X) - 1.00) 

C Compute correlation BT/ M AND GS 
IF (EM. LE .0.10) THEN 
Al = Al 
A2 = A2 
GO TO 160 
ENDIF 

IF (EM. LE .0.65) THEN 

Al = A3 

A2 = A4 

GO TO 160 

ENDIF 

Al = A5 

A2 = A6 

160 CORR = Al - EM*A2 
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nonooo o oonoooooo 


C ESTIMATE CORRELATION BY ITERATION 
EMIN - 1000 
IF (CORR.GE. 0.0) THEN 
CMGS =0.0 
ELSE 

CMGS = -.03 
ENDIF 

Option to compute CMGS by empirical formula (not used) 

IF(EM.GT.0. 48.AND.EM.LT. 0.81) THEN 
E = EM*100 

CMGS = 51.34 - 1 . 957*E + 0 . 0176* (E**2) + 0 . 00003092* (E**3) 
CMGS = CMGS/10000 
GO TO 355 
ENDIF 

Solution by solving quadratic expression 

DO 350 IC = 1,300 

CMGS = CMGS + 0.0001 
ELEFT = (CMGS**2) 

& - (VGI-VM) * (2*VM* (CORR**2) ) /2 

ERIGHT = -CMGS* (2*VM* (CORR**2) ) 

DIFF = ABS (ELEFT - ERIGHT) 

EMIN = MIN (DIFF, EMIN) 

IF (EMIN. EQ. DIFF) THEN 
CM_ACT = CMGS 
ENDIF 
350 CONTINUE 

CMGS = CM_ACT 
355 CONTINUE 

Compute minimization bt/ theoretical and sample reflecctance moments 
DO 100 M=l, 2 

RV (M) = (SER(M) - EGI*ERGI (M) ) /EM 
IF (RV (M) . LT . 0 . 0 ) THEN 
RV (M) =0.00 
ENDIF 

VR (M) = (RV (M) **2) *VM + (EGI**2) *VRGI (M) 

& + (ERGI (M) **2) *VGI + VGI*VRGI(M) 

& - 2*RV(M) *ERGI (M) * (VM + CMGS) 

100 CONTINUE 

CR1R2 = RV (1) *RV (2) *VM 

& - (RV (1) *ERGI (2) + RV (2) *ERGI (1) ) * (VM + CMGS) 

& + (EGI**2) *VRGI (1) + VGI* (VRGI (1) + ERGI (1) *ERGI (2) ) 

NN = NN + 1 
IF (NN. EQ. 10 ) THEN 

WRITE (5, 420) RV(1) , RV (2) , VR(1) ,VR(2) ,CR1R2 
WRITE (5, 420) (Y ( J) , J=l, 5) , VM,F 
WRITE (5, 410) ' RV1 ' , 1 RV2 ' , ’EM', 'AT', 'CMGS', 

& RV ( 1 ) , RV ( 2 ) , EM, AT, CMGS 

410 FORMAT (7X, A, 7X, A, 8X, A, 8X, A, 6X, A, /, 4F10 . 3,F10 . 6) 

420 FORMAT (7E10. 2) 
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c 

c 


c 


c 


NN - 0 
ENDIF 

Y (D “ (VR (1 ) / SVR (1 ) - 1) **2 
Y( 2) = (VR(2) /SVR (2) - 1)**2 
Y (3) = (CR1R2/SCR1R2 - 1)**2 
F = (W1*Y(D) + (W2*Y (2) ) + 

WRITE (5, 420) (Y ( J) , J=1 , 3) , F 


(W3*Y (3) ) 


RETURN 

END 


SUBROUTINE READDAT (AP, M_SAM, CMGS_SAM, 

^ ERGI , VRGI , SER, SVR Qppipox 

REAL ^^(2I.VRGI(2»,SER(2),STO™;fs^, 
& CR1R2 , M_SAM, CMGS SAM 

INTEGER IND ~ 


IF(IND.GT.O) GO TO 100 
C Read fixed parameters 


READ(11,*) AP 

READ (10,-) ERGI (1 ) , ERGI (2), VRGI(l), VRGI (2) 


C Read sample moments 

100 READ (10,*) M_SAM, CMGS_SAM, SER ( 1 ) , SER ( 2 ) 

IND = 1 


SVR(l) ,SVR(2) , SCR1R2 


RETURN 

END 
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INVERSE PROCEDURE FOR CASE V, METHOD 2 

a.) ESTIMATE SOIL BACKGROUND COVER, GI 
(VB1 G_V,FOR) 

Program to compute mean,st. dev., min, and max of a DRAW output 
file for a specified line of R1 and R2 and also compute 
estimate of percent ILLUMINATED GROUND cover, PGIV_EST, for Case V 
simulations and compare it to mean of actual percent cover, PGIV_MEAN. 

REAL RNODE (2) , PGIV_ACT, BIGIR, SMALLIR, BIGVS, SMALLVS, 

& SMALLPGIV, BIGPGIV, 

& RIR_SUM, MEANIR, SDEV_VSSL, SDEV_IRSL, 

& RIRSQ_SUM, RVS_SUM, MEANVS , RVSSQ_SUM, 

& BANDWIDTH,A,B,CVPARM(10) , 

& PGIV_EST_VS , PGIV_EST_IR, PARM, 

& PGIV_MEAN, PGIV_SDEV, PGIV_SUM, PGIV_SSQ, NUM, 

& MACT(IOO) ,MESTVS,MESTIR 

INTEGER NUMPAIRS, NN 
OPEN (10, FILE= ’ G_V .IN', STATUS= ' OLD ' ) 

C ASSIGN INPUT FILE TO FOROll 

OPEN (12, FILE= ' G_V . OUT ' , STATUS 3 'NEW' ) 

OPEN ( 1 3 , FILE= ' G_VL . OUT ’ , STATUS 3 1 NEW * ) 

SMALLIR 3 500.0 
SMALLVS 3 500.0 
SMALLPGIV = 500.0 
NUM 3 0.0 
RIR_SUM 3 0.0 
RIRSQ_SUM = 0.0 
RVS_SUM 3 0.0 
RVSSQ_SUM 3 0.0 
PGIV_SSQ = 0.0 
PGIV_SUM 3 0.0 
NN = 0 

C Read input parameters 

READ (10,*) NUMPAIRS, A, B, BANDWIDTH, PARM 
C Read sdev of soil line and ETA 

READ ( 10 , * ) SDEV_VSSL, SDEV_IRSL 
READ (10,*) ETA 
C Skip header 

READ (11, 300) I PLOT 
300 FORMAT (/, 14, ////////////) 

C Iterate for each data pair 
DO 100 I 3 1, NUMPAIRS 

READ (11,*) (RNODE (K) , K=l, 2) , (CVPARM(J) , J=l, 10) 

PGIV_ACT 3 CVP ARM (PARM) 

C Compute perpendicular distance between a chosen point and 
C the line RIR 3 A*RVIS + B 

AP 3 TAN (-(3. 142/2 - ATAN (A) ) ) 

BP 3 RNODE (2) - (AP*RNODE (1) ) 

RVS 3 (B - BP) / (AP - A) 

RIR 3 A*RVS + B 

DIS =SQRT ( (RIR- RNODE (2) ) **2 + (RVS - RNODE (1 )) **2) 

C WRITE (5,*) AP, BP, RIR, RVS, DIS 

IF (DIS. GT. BANDWIDTH) GO TO 100 
NUM 3 NUM + 1.0 
NN 3 NN + 1 
MACT (NN) 3 CVPARM(l) 


C 

C 

C 

C 

C 

c 

c 

c 

c 
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C Sum IR data 


RIR_SUM = RIR_SUM + RNODE (2) 

RIRSQ_SUM =■ RIRSQ_SUM + (RNODE (2) **2) 

BIGIR = MAX (BIGIR, RNODE (2) ) 

SMALLIR = MIN (SMALLIR, RNODE (2) ) 

C Sum Visible data 

RVS_SUM = RVS_SUM + RNODE (1) 

RVSSQ_SUM = RVSSQ_SUM + (RNODE (1) **2) 

BIGVS = MAX (BIGVS, RNODE (1) ) 

SMALLVS = MIN ( SMALL VS , RNODE ( 1 ) ) 

C Sum percent canopy cover 

PGIV_SUM = PGIV_SUM + PGIV ACT 
PGrv_SSQ = PGIV_SSQ + PGIV - ACT**2 
BIGPGIV = MAX (BIGPGIV, PGIV ACT) 

100 ■ raN <SM«iP<lIV,PGIV_ACTI 

MEANIR = RIR_SUM/NUM 

S : ' N ™*(MEANIR.. 2 ),/ (mjM - 1)1 

C Wri1-I G ?\^f^ EV Z SQRT( ( p GIV SSQ - NUM* (PGIV MEAN**2) ) / (NUM-1) ) 
C Write input and output parameters ~ in 

TO^TE <5 2 210?’ ’S' ;«? ' ; NUM, a, B, BANDWIDTH 

210 r r & 7 xT™ , t ix 6 ;: ;s: ” th 

C Write summary of reflectance statistics 

WRITE (12, 200) 'BAND', 'MEAN', 'SDEV', 'SD/MN', 'MAX' 'MIN' 

& SDEVIR ' SDEVIR/MEANIR, BIGIR, SMALLIR, 

WRITE ( S S - SDEWS ' SDEWS/MEANVS, BIGVS, SMALLVS 

WRITE (5, 200) BAND', 'MEAN', 'SDEV', 'SD/MN', 'MAX’, 'MIN' 

t SDEVIR/MEANIR, BIGIR, SMALLIR, 

200 FORMAT (/, 4X,f, 5^5X^f ' 

„ „ & 6X,A,5F8.2,/,5X,A, 5F8.2) 

- Compute canopy cover using variances 

: COMPUTE E{GI} 

XGIVS = (SDEWS/ SDEV VSSL) 

XGIIR = (SDEVIR/SDEV IRSL) 

PGIV_EST_VS = 100* (XGIVS) 

PGIV_EST_IR = 100* (XGIIR) 

- Compute M based on assumed value of ETA 

MESTVS = 100* (1 - (XGIVS** (1/ (ETA+1) ) ) ) 

. r . MESTIR - 100* (1 - (XGIIR** (1/ (ETA+1) )) ) 

• Write summary of canopy cover statistics 

WRITE (5, 301) 'GI_ACT', 'SDEV', 'G_EST( VIS) ', 'G EST(IR)' 

WRITE (12,301) 

301 

DO 130 J=1,NN 

... WRITE (13, 430) J,MACT(J) , MESTVS, MESTIR 
430 FORMAT (18, 4F8. 2) 

130 CONTINUE 
STOP 
END 
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INVERSE PROCEDURE FOR CASE V, METHOD 2 

b.) ESTIMATE ETA 
(ETA. FOR) 

Program to estimate ETA for large S by selecting 
four conditional lines. 

Use input from M_II.FOR to obtain R(I) and GI(I), RSOIL is 
average reflectance of Soil Line. 

REAL ETA, ETA_INIT, M ( 4 ) , GI (4) , GS (4 ) , R (4) , RSOIL, 

& LHS , RHS 

OPEN (12, FILE= ' ETA .IN’, STATUS= ' OLD ' ) 

OPEN (13, FILE= ' ETA . OUT ' , STATUS= ' NEW ' ) 

READ (12,*) RSOIL 
READ (12, *) (R(I), 1-1,4) 

READ (12, *) (GI(I) ,1=1,4) 

READ (12,*) NETA, ETA_INIT , DEL_ETA 
ETA = ETA_INIT 

WRITE (5, 400) 'A', 'B', 1 C ' , 'D', 'ETA' , ’A-C' , ’B-D' , ' LHS-RHS ' 

WRITE (13, 400) ' A' , 'B', 'C', ' D ' , 'ETA' , 'A-C' , 'B-D' , 'LHS-RHS' 
400 FORMAT (4 (6X, A) , 3 (6X, A) , 2X, A, /) 

C Iterate over ETA 

DO 100 J=1 , NETA 

ETA = ETA + DEL_ETA 

C Compute GS and M 
DO 110 1=1,4 

M(I) = 1.0 - ( (GI (I) ) ** (1/ (ETA+1 . 0) ) ) 

GS (I) = 1.0 - M(I) - GI (I) 

110 CONTINUE 
C Compute LHS 

A = ( (R(l) /M(l) ) - (R(2)/M(2) ))/ 

& ( (GS(1)/M(1) ) - (GS (2) /M(2) ) ) 

B = ( (R(3) /M(3) ) - (R ( 4 ) /M(4) ) ) / 

& ( (GS (3) /M(3) ) - (GS (4) /M(4) ) ) 

LHS = A - B 
C Compute RHS 

C — ( ( (GI(1)/M{1) ) - (GI (2) /M(2) ) ) / 

& ((GS(1)/M(1)) - (GS(2)/M(2) ) ) )* RSOIL 

D = ( ( (GI (3) /M(3) ) - (GI (4)/M(4) ) )/ 

& ( (GS(3) /M(3) ) - (GS(4)/M(4) ) ) )*RSOIL 

RHS = C - D 

DIFF1 = A-C 
DIFF2 = B-D 
QUO = LHS /RHS 

WRITE (5, 410) A, B, C, D, ETA, DIFF1, DIFF2, (LHS-RHS) 

WRITE (13, 410) A, B, C, D, ETA, DIFF1, DIFF2, (LHS-RHS) 

100 CONTINUE 

410 FORMAT (4 (X, F6 . 1) , 4 (X, F8 . 4) ) 

STOP 

END 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 
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c 

c 

c 

c 

c 

c 

c 


INVERSE PROCEDURE FOR BEAVER CREEK 
a.) ESTIMATE BULK PARAMETERS FOR ENTIRE SCENE 
(M_BCA. FOR) 


Estimates 
8 unknowns 


parameters of entire data 
' M, GI, GS, RVRED, RVIR, RSRED 


set from moment analysis 
, RSIR, RGRED, RGIR 


REA1 EVIR ' ^d.rgib, scored, sdgir, 

MEANIR, maxir, minir, rirsq sum 
RRED SUM, MEANRED,MAXRED,MINRED RREDSO qttm 
BANDWIDTH, A, B , VARRED , VAMR^fI, S ‘ 


INTEGER N,NUM, 11,12, NREF,NDAT 


c 


c 


ASSIGN FOR010.DAT BC623B2.0UT 

FILE= ' M _BC. IN r , STATUS = 1 OLD 1 ) 
OPEN (12, FILE= r M BC . OUT 1 , STATUS * 1 NEW 1 ) 
OPEN (1 3, FILE= ■ SLJBC . OUT 1 , STATUS= *NEW f , 


C 

c 

442 


NS = 7 

WRITE (13, 442) ’M_BC, SOIL LINE* , NS, 'NUM' , 'TIME* 

„ 1 f ' ORN 1 , * RED 1 , 1 NIR 1 1 TMP 9 ' 

FORMAT (2X,A,/ ,14, 7 ( / , 2X, A) ) ' ' ™* 

MINIR =500 
MINRED = 500 
NUM = 0 


RIR_SUM = 0 
RIRSQ_SUM = 0 
RRED_SUM = 0 
RREDSQ SUM = 0 


C 


and SOil line Peters 

HEAD (11,*) NDAT 
HEAD (11,*) 11,12 
HEAD (11,*) ALPHA, GAMMA 
HEAD (11,*) FI, AETA 


READ (11,*) RGRED, RGIR, SDGRED, SDGIR 
READ (11,*) A, B, BANDWIDTH 


C Skip header 
READ (10,*) 

HEAD (10,*) ICOL 
DO 100 1=1, ICOL 
READ (10,*) 

100 CONTINUE 


C Read reflectance data 
DO 110 1=1, NDAT 


READ (10, 


r ) N, TIME, (REF (J) ,J=1,NREF) 


c SmS U J 6 conditional moments of arbitrary line A*RR£d + r 
I sISt ^a e pLrs lati0 " bet “ ee " RIR and RVS ' * A-RVIS * B 
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noon 


C New Version 

AP = TAN (- (3 . 142/2-ATAN (A) ) ) 

BP - REF (12) - (AP*REF (II) ) 

RVS = (B-BP) / (AP-A) 

RIR * A* RVS + B 

DIS = SQRT ( (RIR-REF (12) ) **2 + (RVS-REF (II) ) **2) 

IF (DIS.GT. BANDWIDTH) GO TO 110 
NUM = NUM + 1 

C WRITE (13, 441 ) NUM, TIME, (REF (K) , K=l, NREF) 

441 FORMAT (3X, 15, 2X, F10 . 4, 5 (X, F6 . 2) ) 

C Sum IR data 

RIR_SUM = RIR_SUM + REF (12) 

RIRSQ_SUM = RIRSQ_SUM + (REF (12) **2) 

MAXIR = MAX (MAXIR, REF (12) ) 

MINIR = MIN (MINIR, REF (12) ) 

C Sum Visible data 

RRED_SUM = RRED_SUM + REF (II) 

RREDSQ_SUM = RREDSQ_SUM + (REF (II) +*2) 

MAXRED = MAX (MAXRED , REF (II)) 

MINRED = MIN (MINRED, REF (II) ) 

110 CONTINUE 

MEANIR = RIR_SUM/NUM 

SDEVIR = SQRT (RIRSQ_SUM/NUM - (MEANIR** 2) ) 

MEANRED = RRED_SUM/NUM 

SDEVRED = SQRT (RREDSQ_SUM/NUM - (MEANRED** 2) ) 

C Write input and output of conditional line parameters 

WRITE (12,210) 'NUM' , 'A' , 'B' , 'BANDWIDTH' , 'E [GI_RED] ' , ’E[GI_IR] ' , 
& ' SD [GI_RED] ', ' SD [GI_IR] ', NUM, A, B, BANDWIDTH, RGRED, RGIR, 

& SDGRED, SDGIR 

WRITE (5,210) 'NUM' , 'A' , 'B' , 'BANDWIDTH' , ’ E [GI_RED] ’ , 'E(GI_IR] ' , 

& ' SD [GI_RED] ', ' SD [GI_IR] ', NUM, A, B, BANDWIDTH, RGRED, RGIR, 

& SDGRED, SDGIR 

210 FORMAT ( 5X, A, 7X, A, 7X, A, 2X, A,X, A, 2X, A, X, A, X, A, /, 

& 18, 2F8 . 2, F10 . 3, 4F10 . 3) 

WRITE (12, 200) 'BAND', 'MEAN', 'SDEV', 'SD/MN', 'MAX' , 'MIN' , 'ETA' , 

& ' IR ' , MEANIR, SDEVIR, SDEVIR/MEANIR, MAXIR, MINIR, AETA, 

& 'RED' , MEANRED, SDEVRED, SDEVRED /MEANRED, MAXRED , MINRED 

WRITE (5, 200) 'BAND', 'MEAN', 'SDEV', 'SD/MN', 'MAX', 'MIN', 'ETA', 

& ' IR' , MEANIR, SDEVIR, SDEVIR/MEANIR, MAXIR, MINIR, AETA, 

& 'RED' , MEANRED, SDEVRED, SDEVRED /MEANRED, MAXRED, MINRED 

200 FORMAT (/, 4X, A, 4X, A, 4X, A, 3X, A, 5X, A, 5X, A, 5X, A, //, 

& 6X,A, 6F8.3,/,5X,A, 5F8.3) 


Compute parameters of scene 


C Call minimization routine 

CALL EST ( IB , AETA, FI , RGRED , RGIR, SDGRED , SDGIR, MEANRED , MEANIR, 

& SDEVRED, SDEVIR, M, GI,GS, 

& RVRED, RVIR, RSRED, RSIR) 

C Write estimate values 

WRITE (5, 450) 'M', 'GI','GS', ' RV_RED ' , 'RV_IR', ' RS_RED ' , 'RS_IR', 

& M,GI,GS, RVRED, RVIR, RSRED, RSIR 

WRITE (12, 450) 'M', 'GI', 'GS', ' RV_RED ' , ’RV_IR', ' RS_RED ' , 'RS_IR', 

& M,GI,GS, RVRED, RVIR, RSRED, RSIR 

450 FORMAT (/, 9X,A, 8X, A, 8X,A, 4X, A, 5X, A, 4X, A, 5X,A, /, 3 (F10 . 3) , 
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2 (F10 . 3) , 2F10 . 3) 

STOP 

END 


C — 


& 

& 


SUBROUTINE EST (IB, AETA, FI, RGRED, RGIR, SDGRED 
SDEVRED, SDEVIR,M, GI, GS, 

RVRED, RVIR, RSRED, RSIR) 


SDGIR, MEANRED , MEANIR 


& 

Sc 


REAL AETA, FI , RGRED , RGIR, SDGRED , SDGIR, 
SDEVIR,M,GI,GS, 

RVRED, RVIR, RSRED, RSIR 


MEANRED, MEANIR, SDEVRED, 


500 


GI = 0.00 
DIFMIN=100000 . 

DO 500 1=1, 99 
GI = GI + 0.01 

= 1.0/ (1.0 + AETA) 

M = 1.0 - GI** (EXP) 

GS = 1 . 0 - GI - M 

RVRED = (MEANRED - (GI*RGRED) ) / (M + (F1*GS) ) 

RVIR. (MEANIR - (GI*RGIR) ) / (M + (F1*GS) ) ' 

VARRED = SDEVRED**2 1 ° S) ’ 

VARIR = SDEVIR**2 
VARGRED = SDGRED* *2 
VARGIR = SDGIR**2 

VAR_IR_EST = (GI**2) * VARGIR + ( ( (RVRED-RGRED) / (RVIR-RGIR) ) **2) 
* ( (VARRED**2) - (GI**2)*VARGRFn n 2) 

DIFF = ABS (VARIR - VAR IR EST) VARGRED > 

DIFMIN = MIN (DIFF, DIFMIN) ~ 

IF (DIFF . EQ . DIFMIN) THEN 
GI_EST = GI 
EM_EST = M 
GS_EST = GS 
RVRED_EST = RVRED 
RVIRJEST = RVIR 
ENDIF 
CONTINUE 


GI = GI_EST 
M = EM_EST 
GS = GS_EST 
RVRED = RVRED_EST 
RVIR = RVIR EST 
RETURN ~ 

END 
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INVERSE PROCEDURE FOR BEAVER CREEK 
b . ) ESTIMATE BULK PARAMETERS FOR EACH PIXEL 
(M_3C2 . FOR) 

Estimates M directly for each pixel in a conditional parallel line 
Assumes Rgs =0.0 

REAL M, GI , RVRED , RVIR, RGRED , RGIR, 

& TIME, REF (5), 

& BANDWIDTH, A, B, VARRED, VARIR, FI, AETA 

INTEGER N,NUM, 11,12, NREF,NDAT 

C ASSIGN FOR010.DAT BC623B3A.OUT 

OPEN (11, FILE=' M_BC2. IN' , STATUS ='OLD') 

OPEN (12, F ILE='M_BC2. OUT' , STATUS = 'NEW') 

NS = 7 
NUM = 0 

C Read fixed parameters and soil line parameters 
READ (11, *) NREF 
READ (11, *) NDAT 
READ (11, *) 11,12 
READ (11, *) ALPHA, GAMMA 
READ (11,*) FI, AETA 
READ (11,*) RGRED, RGIR, SDGRED , SDGIR 
READ (11,*) A, B, BANDWIDTH 
READ (11,*) RVRED, RVIR 


C Skip header 
READ (10,*) 

READ (10,*) ICOL 
DO 100 1=1, ICOL 
READ(10,*) 

100 CONTINUE 

NNCOL = 9 

WRITE ( 5 , * ) ' M_BC2 ESTIMATES ' 

WRITE (5,*) NNCOL 

WRITE (5, 450) 'N', 'TIME', 'RED', ' IR’ , ' E [M] ' , ' E [GI] ' , ' E [GS] ' , 

& ' Rg [RED ] ' , 1 Rg [ IR] ' 

WRITE (12,*) ' M_BC2 ESTIMATES' 

WRITE (12,*) NNCOL 

WRITE (12, 450) 'N', 'TIME', 'RED', ' IR' , 'E [M] ’ , ' E [GI] ’ , ' E [GS] ' , 

& 'Rg[RED] ' , 'Rg[IR] ' 

450 FORMAT (7X, A, /, 4X,A, /, 5X, A, /, 6X,A, /, 4X,A, /, 3X, A, / , 3X, 

& A, /,X, A, /, 2X,A) 

C Read reflectance data 

DO 110 1=1, NDAT 

READ(10,*) N,TIME, (REF (J),J=1, NREF) 

RED = REF (3) 

RIR = REF (4) 

EM = 0.0 
DIF_MIN=500 
DO 120 J=l,100 
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EM = EM + 0.01 

RIR — ' rEST = (ALPHA*RED) + ( (RVIR-RVRED) *EM) 
+ GAMMA* ( (1-EM) ** (AETA+1 ) ) 

DIFF = ABS (RIR_TEST-RIR) 

DIF_MIN = MIN (DIFF, DIF_MIN) 

IF (DIFF . EQ.DIF_MIN) THEN 
RIR_EST = RIR_TEST 
RED_EST = RED 
EM_EST = EM 

GI_EST = (1-EM) ** (AETA+1) 

GS_EST = 1 - EM - GI_EST 
RG_IR_EST = (RIR - (RVIR+EM) ) /GI EST 

ENDIF^* 0- = (RE ° ' (RVRED * EM > > /GI_EST 
120 CONTINUE 
C Write estimate values 

‘ WRITE<5 ' 460 ’ 

uo‘ ™' 460> 

FORMAT (18, F8 .4,7 (F8 . 3) ) 


460 


STOP 

END 
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